LCOV - code coverage report
Current view: top level - src - qs_wf_history_methods.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:591cf04) Lines: 91.7 % 1281 1175
Test Date: 2026-09-21 02:17:57 Functions: 95.0 % 20 19

            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 Storage of past states of the qs_env.
      10              : !>      Methods to interpolate (or actually normally extrapolate) the
      11              : !>      new guess for density and wavefunctions.
      12              : !> \note
      13              : !>      Most of the last snapshot should actually be in qs_env, but taking
      14              : !>      advantage of it would make the programming much convoluted
      15              : !> \par History
      16              : !>      02.2003 created [fawzi]
      17              : !>      11.2003 Joost VandeVondele : Implemented Nth order PS extrapolation
      18              : !>      02.2005 modified for KG_GPW [MI]
      19              : !> \author fawzi
      20              : ! **************************************************************************************************
      21              : MODULE qs_wf_history_methods
      22              :    USE bibliography,                    ONLY: Kolafa2004,&
      23              :                                               Kuhne2007,&
      24              :                                               VandeVondele2005a,&
      25              :                                               cite_reference
      26              :    USE cell_types,                      ONLY: cell_type,&
      27              :                                               pbc,&
      28              :                                               real_to_scaled
      29              :    USE cp_cfm_basic_linalg,             ONLY: cp_cfm_column_scale,&
      30              :                                               cp_cfm_gemm,&
      31              :                                               cp_cfm_scale_and_add,&
      32              :                                               cp_cfm_trace,&
      33              :                                               cp_cfm_triangular_multiply
      34              :    USE cp_cfm_cholesky,                 ONLY: cp_cfm_cholesky_decompose
      35              :    USE cp_cfm_diag,                     ONLY: cp_cfm_heevd
      36              :    USE cp_cfm_types,                    ONLY: &
      37              :         cp_cfm_create, cp_cfm_get_info, cp_cfm_get_submatrix, cp_cfm_release, &
      38              :         cp_cfm_set_submatrix, cp_cfm_to_cfm, cp_cfm_to_fm, cp_cfm_type, cp_fm_to_cfm
      39              :    USE cp_control_types,                ONLY: dft_control_type
      40              :    USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
      41              :                                               dbcsr_copy,&
      42              :                                               dbcsr_deallocate_matrix,&
      43              :                                               dbcsr_get_info,&
      44              :                                               dbcsr_multiply,&
      45              :                                               dbcsr_p_type,&
      46              :                                               dbcsr_release,&
      47              :                                               dbcsr_type
      48              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_frobenius_norm,&
      49              :                                               dbcsr_trace
      50              :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply,&
      51              :                                               dbcsr_allocate_matrix_set,&
      52              :                                               dbcsr_deallocate_matrix_set
      53              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_scale,&
      54              :                                               cp_fm_scale_and_add
      55              :    USE cp_fm_pool_types,                ONLY: cp_fm_pool_p_type,&
      56              :                                               cp_fm_pool_type,&
      57              :                                               fm_pool_get_el_struct,&
      58              :                                               fm_pools_create_fm_vect,&
      59              :                                               fm_pools_give_back_fm_vect
      60              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      61              :                                               cp_fm_struct_equivalent,&
      62              :                                               cp_fm_struct_release,&
      63              :                                               cp_fm_struct_type
      64              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      65              :                                               cp_fm_get_info,&
      66              :                                               cp_fm_get_submatrix,&
      67              :                                               cp_fm_release,&
      68              :                                               cp_fm_set_all,&
      69              :                                               cp_fm_set_submatrix,&
      70              :                                               cp_fm_to_fm,&
      71              :                                               cp_fm_type
      72              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      73              :                                               cp_logger_type,&
      74              :                                               cp_to_string
      75              :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      76              :                                               cp_print_key_unit_nr,&
      77              :                                               low_print_level
      78              :    USE input_constants,                 ONLY: &
      79              :         wfi_aspc_nr, wfi_frozen_method_nr, wfi_gext_proj_nr, wfi_gext_proj_qtr_nr, &
      80              :         wfi_linear_p_method_nr, wfi_linear_ps_method_nr, wfi_linear_wf_method_nr, &
      81              :         wfi_ps_method_nr, wfi_use_guess_method_nr, wfi_use_prev_p_method_nr, &
      82              :         wfi_use_prev_wf_method_nr
      83              :    USE kinds,                           ONLY: dp
      84              :    USE kpoint_types,                    ONLY: get_kpoint_info,&
      85              :                                               kpoint_env_type,&
      86              :                                               kpoint_type
      87              :    USE mathconstants,                   ONLY: twopi,&
      88              :                                               z_one,&
      89              :                                               z_zero
      90              :    USE mathlib,                         ONLY: binomial
      91              :    USE message_passing,                 ONLY: mp_para_env_type
      92              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      93              :    USE particle_types,                  ONLY: particle_type
      94              :    USE pw_env_types,                    ONLY: pw_env_get,&
      95              :                                               pw_env_type
      96              :    USE pw_methods,                      ONLY: pw_copy
      97              :    USE pw_pool_types,                   ONLY: pw_pool_type
      98              :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      99              :                                               pw_r3d_rs_type
     100              :    USE qs_density_matrices,             ONLY: calculate_density_matrix
     101              :    USE qs_environment_types,            ONLY: get_qs_env,&
     102              :                                               qs_environment_type,&
     103              :                                               set_qs_env
     104              :    USE qs_kpoint_operators,             ONLY: kpoint_operator_get_all,&
     105              :                                               kpoint_slot_s,&
     106              :                                               kpoint_spin_free
     107              :    USE qs_kpoint_state,                 ONLY: qs_kpoint_state_commit
     108              :    USE qs_ks_types,                     ONLY: qs_ks_did_change
     109              :    USE qs_matrix_pools,                 ONLY: mpools_get,&
     110              :                                               qs_matrix_pools_type
     111              :    USE qs_mo_methods,                   ONLY: make_basis_cholesky,&
     112              :                                               make_basis_lowdin,&
     113              :                                               make_basis_simple,&
     114              :                                               make_basis_sm
     115              :    USE qs_mo_types,                     ONLY: get_mo_set,&
     116              :                                               mo_set_type
     117              :    USE qs_rho_methods,                  ONLY: qs_rho_update_rho
     118              :    USE qs_rho_types,                    ONLY: qs_rho_get,&
     119              :                                               qs_rho_type
     120              :    USE qs_scf_types,                    ONLY: ot_method_nr,&
     121              :                                               qs_scf_env_type
     122              :    USE qs_wf_history_types,             ONLY: qs_wf_history_type,&
     123              :                                               qs_wf_snapshot_type,&
     124              :                                               wfi_get_snapshot,&
     125              :                                               wfi_release
     126              :    USE scf_control_types,               ONLY: scf_control_type
     127              : #include "./base/base_uses.f90"
     128              : 
     129              :    IMPLICIT NONE
     130              :    PRIVATE
     131              : 
     132              :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
     133              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_wf_history_methods'
     134              : 
     135              :    PUBLIC :: wfi_create, wfi_update, wfi_create_for_kp, &
     136              :              wfi_extrapolate, wfi_get_method_label, &
     137              :              reorthogonalize_vectors, wfi_purge_history
     138              : 
     139              : CONTAINS
     140              : 
     141              : ! **************************************************************************************************
     142              : !> \brief allocates and initialize a wavefunction snapshot
     143              : !> \param snapshot the snapshot to create
     144              : !> \par History
     145              : !>      02.2003 created [fawzi]
     146              : !>      02.2005 added wf_mol [MI]
     147              : !> \author fawzi
     148              : ! **************************************************************************************************
     149        12176 :    SUBROUTINE wfs_create(snapshot)
     150              :       TYPE(qs_wf_snapshot_type), INTENT(OUT)             :: snapshot
     151              : 
     152              :       NULLIFY (snapshot%wf, snapshot%rho_r, &
     153              :                snapshot%rho_g, snapshot%rho_ao, snapshot%rho_ao_kp, &
     154              :                snapshot%overlap, snapshot%wf_kp, snapshot%overlap_cfm_kp, &
     155              :                snapshot%kp_pbc_shift, snapshot%rho_frozen)
     156        12176 :       snapshot%dt = 1.0_dp
     157        12176 :    END SUBROUTINE wfs_create
     158              : 
     159              : ! **************************************************************************************************
     160              : !> \brief updates the given snapshot
     161              : !> \param snapshot the snapshot to be updated
     162              : !> \param wf_history the history
     163              : !> \param qs_env the qs_env that should be snapshotted
     164              : !> \param dt the time of the snapshot (wrt. to the previous snapshot)
     165              : !> \par History
     166              : !>      02.2003 created [fawzi]
     167              : !>      02.2005 added kg_fm_mol_set for KG_GPW [MI]
     168              : !> \author fawzi
     169              : ! **************************************************************************************************
     170        22616 :    SUBROUTINE wfs_update(snapshot, wf_history, qs_env, dt)
     171              :       TYPE(qs_wf_snapshot_type), POINTER                 :: snapshot
     172              :       TYPE(qs_wf_history_type), POINTER                  :: wf_history
     173              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     174              :       REAL(KIND=dp), INTENT(in), OPTIONAL                :: dt
     175              : 
     176              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'wfs_update'
     177              : 
     178              :       INTEGER                                            :: handle, ic, ikp, img, ispin, kplocal, &
     179              :                                                             nc, nimg, nspin_kp, nspins
     180              :       INTEGER, DIMENSION(2)                              :: kp_range
     181              :       TYPE(cell_type), POINTER                           :: cell
     182        22616 :       TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER     :: ao_mo_pools
     183              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     184        22616 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s, rho_ao
     185        22616 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: rho_ao_kp
     186              :       TYPE(dft_control_type), POINTER                    :: dft_control
     187              :       TYPE(kpoint_env_type), POINTER                     :: kp
     188              :       TYPE(kpoint_type), POINTER                         :: kpoints
     189        22616 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     190        22616 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     191        22616 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g
     192              :       TYPE(pw_env_type), POINTER                         :: pw_env
     193              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     194        22616 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r
     195              :       TYPE(qs_rho_type), POINTER                         :: rho
     196              : 
     197        22616 :       CALL timeset(routineN, handle)
     198              : 
     199        22616 :       NULLIFY (pw_env, auxbas_pw_pool, ao_mo_pools, dft_control, mos, mo_coeff, &
     200        22616 :                rho, rho_r, rho_g, rho_ao, matrix_s, kpoints, kp, cell, particle_set)
     201              :       CALL get_qs_env(qs_env, pw_env=pw_env, &
     202        22616 :                       dft_control=dft_control, rho=rho, cell=cell, particle_set=particle_set)
     203        22616 :       CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_pools)
     204        22616 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     205              : 
     206        22616 :       CPASSERT(ASSOCIATED(wf_history))
     207        22616 :       CPASSERT(ASSOCIATED(dft_control))
     208        22616 :       IF (.NOT. ASSOCIATED(snapshot)) THEN
     209        12176 :          ALLOCATE (snapshot)
     210        12176 :          CALL wfs_create(snapshot)
     211              :       END IF
     212        22616 :       CPASSERT(wf_history%ref_count > 0)
     213              : 
     214        22616 :       nspins = dft_control%nspins
     215        22616 :       snapshot%dt = 1.0_dp
     216        22616 :       IF (PRESENT(dt)) snapshot%dt = dt
     217        22616 :       IF (wf_history%store_wf) THEN
     218        19464 :          CALL get_qs_env(qs_env, mos=mos)
     219        19464 :          IF (.NOT. ASSOCIATED(snapshot%wf)) THEN
     220              :             CALL fm_pools_create_fm_vect(ao_mo_pools, snapshot%wf, &
     221         9818 :                                          name="ws_snap-ws")
     222         9818 :             CPASSERT(nspins == SIZE(snapshot%wf))
     223              :          END IF
     224        41294 :          DO ispin = 1, nspins
     225        21830 :             CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
     226        41294 :             CALL cp_fm_to_fm(mo_coeff, snapshot%wf(ispin))
     227              :          END DO
     228              :       ELSE
     229         3152 :          CALL fm_pools_give_back_fm_vect(ao_mo_pools, snapshot%wf)
     230              :       END IF
     231              : 
     232        22616 :       IF (wf_history%store_rho_r) THEN
     233            0 :          CALL qs_rho_get(rho, rho_r=rho_r)
     234            0 :          CPASSERT(ASSOCIATED(rho_r))
     235            0 :          IF (.NOT. ASSOCIATED(snapshot%rho_r)) THEN
     236            0 :             ALLOCATE (snapshot%rho_r(nspins))
     237            0 :             DO ispin = 1, nspins
     238            0 :                CALL auxbas_pw_pool%create_pw(snapshot%rho_r(ispin))
     239              :             END DO
     240              :          END IF
     241            0 :          DO ispin = 1, nspins
     242            0 :             CALL pw_copy(rho_r(ispin), snapshot%rho_r(ispin))
     243              :          END DO
     244        22616 :       ELSE IF (ASSOCIATED(snapshot%rho_r)) THEN
     245            0 :          DO ispin = 1, SIZE(snapshot%rho_r)
     246            0 :             CALL auxbas_pw_pool%give_back_pw(snapshot%rho_r(ispin))
     247              :          END DO
     248            0 :          DEALLOCATE (snapshot%rho_r)
     249              :       END IF
     250              : 
     251        22616 :       IF (wf_history%store_rho_g) THEN
     252            0 :          CALL qs_rho_get(rho, rho_g=rho_g)
     253            0 :          CPASSERT(ASSOCIATED(rho_g))
     254            0 :          IF (.NOT. ASSOCIATED(snapshot%rho_g)) THEN
     255            0 :             ALLOCATE (snapshot%rho_g(nspins))
     256            0 :             DO ispin = 1, nspins
     257            0 :                CALL auxbas_pw_pool%create_pw(snapshot%rho_g(ispin))
     258              :             END DO
     259              :          END IF
     260            0 :          DO ispin = 1, nspins
     261            0 :             CALL pw_copy(rho_g(ispin), snapshot%rho_g(ispin))
     262              :          END DO
     263        22616 :       ELSE IF (ASSOCIATED(snapshot%rho_g)) THEN
     264            0 :          DO ispin = 1, SIZE(snapshot%rho_g)
     265            0 :             CALL auxbas_pw_pool%give_back_pw(snapshot%rho_g(ispin))
     266              :          END DO
     267            0 :          DEALLOCATE (snapshot%rho_g)
     268              :       END IF
     269              : 
     270        22616 :       IF (ASSOCIATED(snapshot%rho_ao)) THEN ! the sparsity might be different
     271              :          ! (future struct:check)
     272          262 :          CALL dbcsr_deallocate_matrix_set(snapshot%rho_ao)
     273              :       END IF
     274        22616 :       IF (wf_history%store_rho_ao) THEN
     275          310 :          CALL qs_rho_get(rho, rho_ao=rho_ao)
     276          310 :          CPASSERT(ASSOCIATED(rho_ao))
     277              : 
     278          310 :          CALL dbcsr_allocate_matrix_set(snapshot%rho_ao, nspins)
     279          770 :          DO ispin = 1, nspins
     280          460 :             ALLOCATE (snapshot%rho_ao(ispin)%matrix)
     281          770 :             CALL dbcsr_copy(snapshot%rho_ao(ispin)%matrix, rho_ao(ispin)%matrix)
     282              :          END DO
     283              :       END IF
     284              : 
     285        22616 :       IF (ASSOCIATED(snapshot%rho_ao_kp)) THEN ! the sparsity might be different
     286              :          ! (future struct:check)
     287          220 :          CALL dbcsr_deallocate_matrix_set(snapshot%rho_ao_kp)
     288              :       END IF
     289        22616 :       IF (wf_history%store_rho_ao_kp) THEN
     290          232 :          CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
     291          232 :          CPASSERT(ASSOCIATED(rho_ao_kp))
     292              : 
     293          232 :          nimg = dft_control%nimages
     294          232 :          CALL dbcsr_allocate_matrix_set(snapshot%rho_ao_kp, nspins, nimg)
     295          554 :          DO ispin = 1, nspins
     296        34092 :             DO img = 1, nimg
     297        33538 :                ALLOCATE (snapshot%rho_ao_kp(ispin, img)%matrix)
     298              :                CALL dbcsr_copy(snapshot%rho_ao_kp(ispin, img)%matrix, &
     299        33860 :                                rho_ao_kp(ispin, img)%matrix)
     300              :             END DO
     301              :          END DO
     302              :       END IF
     303              : 
     304        22616 :       IF (ASSOCIATED(snapshot%overlap)) THEN ! the sparsity might be different
     305              :          ! (future struct:check)
     306         7038 :          CALL dbcsr_deallocate_matrix(snapshot%overlap)
     307              :       END IF
     308        22616 :       IF (wf_history%store_overlap) THEN
     309        15946 :          CALL get_qs_env(qs_env, matrix_s=matrix_s)
     310        15946 :          CPASSERT(ASSOCIATED(matrix_s))
     311        15946 :          CPASSERT(ASSOCIATED(matrix_s(1)%matrix))
     312        15946 :          ALLOCATE (snapshot%overlap)
     313        15946 :          CALL dbcsr_copy(snapshot%overlap, matrix_s(1)%matrix)
     314              :       END IF
     315              : 
     316        22616 :       CALL get_qs_env(qs_env, kpoints=kpoints)
     317        22616 :       IF (ASSOCIATED(kpoints)) THEN
     318        22616 :          IF (ASSOCIATED(kpoints%kp_env)) THEN
     319              :             ! --- k-point WFN snapshot: store complex MO coefficients per local k-point ---
     320         2826 :             IF (wf_history%store_wf_kp) THEN
     321         2594 :                CALL get_kpoint_info(kpoints, kp_range=kp_range)
     322         2594 :                kplocal = kp_range(2) - kp_range(1) + 1
     323         2594 :                nspin_kp = SIZE(kpoints%kp_env(1)%kpoint_env%mos, 2)
     324         2594 :                nc = SIZE(kpoints%kp_env(1)%kpoint_env%mos, 1) ! 2=complex, 1=real
     325              : 
     326         2594 :                CALL wfi_store_kp_pbc_shift(snapshot, cell, particle_set)
     327              : 
     328         2594 :                IF (ASSOCIATED(snapshot%wf_kp)) THEN
     329          842 :                   DO ikp = 1, SIZE(snapshot%wf_kp, 1)
     330         1926 :                      DO ic = 1, SIZE(snapshot%wf_kp, 2)
     331         2710 :                         DO ispin = 1, SIZE(snapshot%wf_kp, 3)
     332         2168 :                            CALL cp_fm_release(snapshot%wf_kp(ikp, ic, ispin))
     333              :                         END DO
     334              :                      END DO
     335              :                   END DO
     336          300 :                   DEALLOCATE (snapshot%wf_kp)
     337              :                END IF
     338              : 
     339        31438 :                ALLOCATE (snapshot%wf_kp(kplocal, nc, nspin_kp))
     340         7141 :                DO ikp = 1, kplocal
     341         4547 :                   kp => kpoints%kp_env(ikp)%kpoint_env
     342        12298 :                   DO ispin = 1, nspin_kp
     343        20018 :                      DO ic = 1, nc
     344        10314 :                         CALL get_mo_set(kp%mos(ic, ispin), mo_coeff=mo_coeff)
     345              :                         CALL cp_fm_create(snapshot%wf_kp(ikp, ic, ispin), &
     346              :                                           mo_coeff%matrix_struct, &
     347        10314 :                                           name="wfkp_snap")
     348        15471 :                         CALL cp_fm_to_fm(mo_coeff, snapshot%wf_kp(ikp, ic, ispin))
     349              :                      END DO
     350              :                   END DO
     351              :                END DO
     352              :             END IF
     353              : 
     354         2826 :             IF (wf_history%store_overlap_kp) CALL wfs_overlap_kp(snapshot, qs_env)
     355              :          END IF
     356              :       END IF
     357              : 
     358              :       IF (wf_history%store_frozen_density) THEN
     359              :          ! do nothing
     360              :          ! CALL deallocate_matrix_set(snapshot%rho_frozen%rho_ao)
     361              :       END IF
     362              : 
     363        22616 :       CALL timestop(handle)
     364              : 
     365        22616 :    END SUBROUTINE wfs_update
     366              : 
     367              : ! **************************************************************************************************
     368              : !> \brief Store S(k) in the current neighbor-list and PBC convention.
     369              : !>        Reuse the unmodified overlap cache from this SCF run when available.
     370              : !>        Otherwise build the snapshot now, never from S(R) at a later geometry.
     371              : !> \param snapshot owns its matrices independently of the SCF cache and operator context
     372              : !> \param qs_env current environment; its SCF overlap cache is invalidated by init_scf_run
     373              : ! **************************************************************************************************
     374         2580 :    SUBROUTINE wfs_overlap_kp(snapshot, qs_env)
     375              :       TYPE(qs_wf_snapshot_type), INTENT(INOUT)           :: snapshot
     376              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     377              : 
     378              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'wfs_overlap_kp'
     379              : 
     380              :       INTEGER                                            :: handle, ikp, kplocal, use_cache
     381         2580 :       TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER     :: ao_ao_fm_pools
     382              :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct
     383         2580 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s
     384              :       TYPE(kpoint_type), POINTER                         :: kpoints
     385              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     386              :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
     387              : 
     388         2580 :       CALL timeset(routineN, handle)
     389         2580 :       CALL get_qs_env(qs_env, kpoints=kpoints, matrix_s_kp=matrix_s, scf_env=scf_env)
     390         2580 :       para_env => kpoints%blacs_env_all%para_env
     391         2580 :       kplocal = SIZE(kpoints%kp_env)
     392         2580 :       CPASSERT(kplocal > 0)
     393         2580 :       CALL mpools_get(kpoints%mpools, ao_ao_fm_pools=ao_ao_fm_pools)
     394         2580 :       matrix_struct => fm_pool_get_el_struct(ao_ao_fm_pools(1)%pool)
     395              : 
     396              :       ! All ranks must choose the same path if any cache is missing or incompatible.
     397         2580 :       use_cache = 0
     398         2580 :       IF (ALLOCATED(kpoints%scf_smat)) THEN
     399         2558 :          IF (SIZE(kpoints%scf_smat) == kplocal) THEN
     400         2558 :             use_cache = 1
     401         7039 :             DO ikp = 1, kplocal
     402         7039 :                IF (.NOT. cp_fm_struct_equivalent(kpoints%scf_smat(ikp)%matrix_struct, matrix_struct)) use_cache = 0
     403              :             END DO
     404              :          END IF
     405              :       END IF
     406         2580 :       CALL para_env%min(use_cache)
     407              : 
     408         2580 :       IF (ASSOCIATED(snapshot%overlap_cfm_kp)) THEN
     409          822 :          DO ikp = 1, SIZE(snapshot%overlap_cfm_kp)
     410          822 :             CALL cp_cfm_release(snapshot%overlap_cfm_kp(ikp))
     411              :          END DO
     412          296 :          DEALLOCATE (snapshot%overlap_cfm_kp)
     413              :       END IF
     414        12255 :       ALLOCATE (snapshot%overlap_cfm_kp(kplocal))
     415              :       ! Matrix allocation and structure refcounts stay outside OpenMP workers.
     416         7095 :       DO ikp = 1, kplocal
     417         7095 :          CALL cp_cfm_create(snapshot%overlap_cfm_kp(ikp), matrix_struct)
     418              :       END DO
     419              : 
     420         2580 :       IF (use_cache == 1) THEN
     421         7039 :          DO ikp = 1, kplocal
     422         7039 :             CALL cp_cfm_to_cfm(kpoints%scf_smat(ikp), snapshot%overlap_cfm_kp(ikp))
     423              :          END DO
     424              :       ELSE
     425              :          CALL kpoint_operator_get_all(kpoints, scf_env%scf_work1(1), matrix_s, kpoint_slot_s, &
     426           22 :                                       snapshot%overlap_cfm_kp, matrix_row=kpoint_spin_free)
     427              :       END IF
     428         2580 :       CALL timestop(handle)
     429              : 
     430         2580 :    END SUBROUTINE wfs_overlap_kp
     431              : 
     432              : ! **************************************************************************************************
     433              : !> \brief ...
     434              : !> \param wf_history ...
     435              : !> \param interpolation_method_nr the tag of the method used for
     436              : !>        the extrapolation of the initial density for the next md step
     437              : !>        (see qs_wf_history_types:wfi_*_method_nr)
     438              : !> \param extrapolation_order ...
     439              : !> \param has_unit_metric ...
     440              : !> \par History
     441              : !>      02.2003 created [fawzi]
     442              : !> \author fawzi
     443              : ! **************************************************************************************************
     444         9132 :    SUBROUTINE wfi_create(wf_history, interpolation_method_nr, extrapolation_order, &
     445              :                          has_unit_metric)
     446              :       TYPE(qs_wf_history_type), POINTER                  :: wf_history
     447              :       INTEGER, INTENT(in)                                :: interpolation_method_nr, &
     448              :                                                             extrapolation_order
     449              :       LOGICAL, INTENT(IN)                                :: has_unit_metric
     450              : 
     451              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'wfi_create'
     452              : 
     453              :       INTEGER                                            :: handle, i
     454              : 
     455         9132 :       CALL timeset(routineN, handle)
     456              : 
     457         9132 :       ALLOCATE (wf_history)
     458         9132 :       wf_history%ref_count = 1
     459         9132 :       wf_history%memory_depth = 0
     460         9132 :       wf_history%snapshot_count = 0
     461         9132 :       wf_history%last_state_index = 1
     462              :       wf_history%store_wf = .FALSE.
     463              :       wf_history%store_rho_r = .FALSE.
     464              :       wf_history%store_rho_g = .FALSE.
     465              :       wf_history%store_rho_ao = .FALSE.
     466              :       wf_history%store_rho_ao_kp = .FALSE.
     467              :       wf_history%store_overlap = .FALSE.
     468              :       wf_history%store_wf_kp = .FALSE.
     469              :       wf_history%store_overlap_kp = .FALSE.
     470              :       wf_history%store_frozen_density = .FALSE.
     471              :       NULLIFY (wf_history%past_states)
     472              : 
     473         9132 :       wf_history%interpolation_method_nr = interpolation_method_nr
     474              : 
     475              :       SELECT CASE (wf_history%interpolation_method_nr)
     476              :       CASE (wfi_use_guess_method_nr)
     477              :          wf_history%memory_depth = 0
     478              :       CASE (wfi_use_prev_wf_method_nr)
     479           64 :          wf_history%memory_depth = 0
     480              :       CASE (wfi_use_prev_p_method_nr)
     481           64 :          wf_history%memory_depth = 1
     482           64 :          wf_history%store_rho_ao = .TRUE.
     483              :       CASE (wfi_linear_wf_method_nr)
     484            4 :          wf_history%memory_depth = 2
     485            4 :          wf_history%store_wf = .TRUE.
     486              :       CASE (wfi_linear_p_method_nr)
     487            6 :          wf_history%memory_depth = 2
     488            6 :          wf_history%store_rho_ao = .TRUE.
     489              :       CASE (wfi_linear_ps_method_nr)
     490            6 :          wf_history%memory_depth = 2
     491            6 :          wf_history%store_wf = .TRUE.
     492            6 :          IF (.NOT. has_unit_metric) wf_history%store_overlap = .TRUE.
     493              :       CASE (wfi_ps_method_nr)
     494          345 :          CALL cite_reference(VandeVondele2005a)
     495          345 :          wf_history%memory_depth = extrapolation_order + 1
     496          345 :          wf_history%store_wf = .TRUE.
     497          345 :          wf_history%store_wf_kp = .TRUE.
     498          345 :          IF (.NOT. has_unit_metric) THEN
     499          341 :             wf_history%store_overlap = .TRUE.
     500          341 :             wf_history%store_overlap_kp = .TRUE.
     501              :          END IF
     502              :       CASE (wfi_frozen_method_nr)
     503            4 :          wf_history%memory_depth = 1
     504            4 :          wf_history%store_frozen_density = .TRUE.
     505              :       CASE (wfi_aspc_nr)
     506         8317 :          wf_history%memory_depth = extrapolation_order + 2
     507         8317 :          wf_history%store_wf = .TRUE.
     508         8317 :          wf_history%store_wf_kp = .TRUE.
     509         8317 :          IF (.NOT. has_unit_metric) THEN
     510         7335 :             wf_history%store_overlap = .TRUE.
     511         7335 :             wf_history%store_overlap_kp = .TRUE.
     512              :          END IF
     513              :       CASE (wfi_gext_proj_nr)
     514           26 :          wf_history%memory_depth = extrapolation_order
     515           26 :          wf_history%store_wf = .TRUE.
     516           26 :          wf_history%store_wf_kp = .TRUE.
     517           26 :          wf_history%store_overlap = .TRUE.
     518           26 :          wf_history%store_overlap_kp = .TRUE.
     519              :       CASE (wfi_gext_proj_qtr_nr)
     520            6 :          wf_history%memory_depth = extrapolation_order
     521            6 :          wf_history%store_wf = .TRUE.
     522            6 :          wf_history%store_wf_kp = .TRUE.
     523            6 :          wf_history%store_overlap = .TRUE.
     524            6 :          wf_history%store_overlap_kp = .TRUE.
     525              :       CASE default
     526              :          CALL cp_abort(__LOCATION__, &
     527              :                        "Unknown interpolation method: "// &
     528         9132 :                        TRIM(ADJUSTL(cp_to_string(interpolation_method_nr))))
     529              :       END SELECT
     530        69897 :       ALLOCATE (wf_history%past_states(wf_history%memory_depth))
     531              : 
     532        51987 :       DO i = 1, SIZE(wf_history%past_states)
     533        51987 :          NULLIFY (wf_history%past_states(i)%snapshot)
     534              :       END DO
     535              : 
     536         9132 :       CALL timestop(handle)
     537         9132 :    END SUBROUTINE wfi_create
     538              : 
     539              : ! **************************************************************************************************
     540              : !> \brief Adapts wf_history storage flags for k-point calculations.
     541              : !>        For ASPC, switches from Gamma WFN storage to k-point WFN storage.
     542              : !>        Other WFN-based methods remain blocked.
     543              : !> \param wf_history ...
     544              : !> \par History
     545              : !>      06.2015 created [jhu]
     546              : !> \author jhu
     547              : ! **************************************************************************************************
     548          680 :    SUBROUTINE wfi_create_for_kp(wf_history)
     549              :       TYPE(qs_wf_history_type), POINTER                  :: wf_history
     550              : 
     551              :       INTEGER                                            :: i
     552              : 
     553          680 :       CPASSERT(ASSOCIATED(wf_history))
     554          680 :       IF (wf_history%store_rho_ao) THEN
     555           10 :          wf_history%store_rho_ao_kp = .TRUE.
     556           10 :          wf_history%store_rho_ao = .FALSE.
     557              :       END IF
     558              :       ! KP-compatible WFN history: store complex k-point MOs in snapshots.
     559              :       ! USE_PREV_WF needs one snapshot as well, since the PBC image convention
     560              :       ! of the saved WFN has to be known before reorthogonalization.
     561          680 :       IF (wf_history%interpolation_method_nr == wfi_use_prev_wf_method_nr) THEN
     562           30 :          wf_history%memory_depth = 1
     563           30 :          wf_history%store_wf_kp = .TRUE.
     564           30 :          wf_history%store_wf = .FALSE.
     565           30 :          wf_history%store_overlap = .FALSE.
     566           30 :          IF (ASSOCIATED(wf_history%past_states)) DEALLOCATE (wf_history%past_states)
     567          120 :          ALLOCATE (wf_history%past_states(wf_history%memory_depth))
     568           60 :          DO i = 1, SIZE(wf_history%past_states)
     569           60 :             NULLIFY (wf_history%past_states(i)%snapshot)
     570              :          END DO
     571          650 :       ELSE IF (wf_history%store_wf_kp) THEN
     572          504 :          wf_history%store_wf = .FALSE.
     573          504 :          wf_history%store_overlap = .FALSE.
     574              :          ! store_wf_kp and store_overlap_kp remain TRUE
     575              :       ELSE
     576              :          ! Linear methods (except LINEAR_P) are still blocked
     577          146 :          IF (wf_history%store_wf .OR. wf_history%store_overlap) THEN
     578              :             CALL cp_abort(__LOCATION__, &
     579            0 :                           "Linear WFN-based extrapolation methods not implemented for k-points.")
     580              :          END IF
     581              :       END IF
     582          680 :       IF (wf_history%store_frozen_density) THEN
     583            0 :          CPABORT("Frozen density initialization method not possible for kpoints.")
     584              :       END IF
     585              : 
     586          680 :    END SUBROUTINE wfi_create_for_kp
     587              : 
     588              : ! **************************************************************************************************
     589              : !> \brief returns a string describing the interpolation method
     590              : !> \param method_nr ...
     591              : !> \return ...
     592              : !> \par History
     593              : !>      02.2003 created [fawzi]
     594              : !> \author fawzi
     595              : ! **************************************************************************************************
     596        13615 :    FUNCTION wfi_get_method_label(method_nr) RESULT(res)
     597              :       INTEGER, INTENT(in)                                :: method_nr
     598              :       CHARACTER(len=30)                                  :: res
     599              : 
     600        13615 :       res = "unknown"
     601        13853 :       SELECT CASE (method_nr)
     602              :       CASE (wfi_use_prev_p_method_nr)
     603          238 :          res = "previous_p"
     604              :       CASE (wfi_use_prev_wf_method_nr)
     605          333 :          res = "previous_wf"
     606              :       CASE (wfi_use_guess_method_nr)
     607         5177 :          res = "initial_guess"
     608              :       CASE (wfi_linear_wf_method_nr)
     609            2 :          res = "mo linear"
     610              :       CASE (wfi_linear_p_method_nr)
     611            3 :          res = "P linear"
     612              :       CASE (wfi_linear_ps_method_nr)
     613            6 :          res = "PS linear"
     614              :       CASE (wfi_ps_method_nr)
     615          188 :          res = "PS Nth order"
     616              :       CASE (wfi_frozen_method_nr)
     617            4 :          res = "frozen density approximation"
     618              :       CASE (wfi_aspc_nr)
     619         7577 :          res = "ASPC"
     620              :       CASE (wfi_gext_proj_nr)
     621           73 :          res = "GEXT_PROJ"
     622              :       CASE (wfi_gext_proj_qtr_nr)
     623           14 :          res = "GEXT_PROJ_QTR"
     624              :       CASE default
     625              :          CALL cp_abort(__LOCATION__, &
     626              :                        "Unknown interpolation method: "// &
     627        13615 :                        TRIM(ADJUSTL(cp_to_string(method_nr))))
     628              :       END SELECT
     629        13615 :    END FUNCTION wfi_get_method_label
     630              : 
     631              : ! **************************************************************************************************
     632              : !> \brief calculates the new starting state for the scf for the next
     633              : !>      wf optimization
     634              : !> \param wf_history the previous history needed to extrapolate
     635              : !> \param qs_env the qs env with the latest result, and that will contain
     636              : !>        the new starting state
     637              : !> \param dt the time at which to extrapolate (wrt. to the last snapshot)
     638              : !> \param extrapolation_method_nr returns the extrapolation method used
     639              : !> \param orthogonal_wf ...
     640              : !> \par History
     641              : !>      02.2003 created [fawzi]
     642              : !>      11.2003 Joost VandeVondele : Implemented Nth order PS extrapolation
     643              : !>      04.2026 Michele Nottoli : Added GEXT_PROJ and GEXT_PROJ_QTR extrapolations
     644              : !> \author fawzi
     645              : ! **************************************************************************************************
     646        26933 :    SUBROUTINE wfi_extrapolate(wf_history, qs_env, dt, extrapolation_method_nr, &
     647              :                               orthogonal_wf)
     648              :       TYPE(qs_wf_history_type), POINTER                  :: wf_history
     649              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     650              :       REAL(KIND=dp), INTENT(IN)                          :: dt
     651              :       INTEGER, INTENT(OUT), OPTIONAL                     :: extrapolation_method_nr
     652              :       LOGICAL, INTENT(OUT), OPTIONAL                     :: orthogonal_wf
     653              : 
     654              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'wfi_extrapolate'
     655              : 
     656              :       INTEGER                                            :: actual_extrapolation_method_nr, handle, &
     657              :                                                             i, img, io_unit, ispin, k, n, nmo, &
     658              :                                                             nvec, print_level
     659              :       LOGICAL                                            :: do_kpoints, my_orthogonal_wf, use_overlap
     660              :       REAL(KIND=dp)                                      :: alpha, t0, t1, t2
     661        26933 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: coeffs
     662        26933 :       TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER     :: ao_mo_fm_pools
     663              :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct, matrix_struct_new
     664              :       TYPE(cp_fm_type)                                   :: csc, fm_tmp
     665              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     666              :       TYPE(cp_logger_type), POINTER                      :: logger
     667        26933 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s, rho_ao, rho_frozen_ao
     668        26933 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: rho_ao_kp
     669        26933 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     670              :       TYPE(qs_rho_type), POINTER                         :: rho
     671              :       TYPE(qs_wf_snapshot_type), POINTER                 :: t0_state, t1_state
     672              : 
     673        26933 :       NULLIFY (mos, ao_mo_fm_pools, t0_state, t1_state, mo_coeff, &
     674        26933 :                rho, rho_ao, rho_frozen_ao)
     675              : 
     676        26933 :       use_overlap = wf_history%store_overlap
     677              : 
     678        26933 :       CALL timeset(routineN, handle)
     679        26933 :       logger => cp_get_default_logger()
     680        26933 :       print_level = logger%iter_info%print_level
     681              :       io_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%SCF%PRINT%PROGRAM_RUN_INFO", &
     682        26933 :                                      extension=".scfLog")
     683              : 
     684        26933 :       CPASSERT(ASSOCIATED(wf_history))
     685        26933 :       CPASSERT(wf_history%ref_count > 0)
     686        26933 :       CPASSERT(ASSOCIATED(qs_env))
     687        26933 :       CALL get_qs_env(qs_env, mos=mos, rho=rho, do_kpoints=do_kpoints)
     688        26933 :       CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
     689              :       ! chooses the method for this extrapolation
     690        26933 :       IF (wf_history%snapshot_count < 1) THEN
     691              :          actual_extrapolation_method_nr = wfi_use_guess_method_nr
     692              :       ELSE
     693        16742 :          actual_extrapolation_method_nr = wf_history%interpolation_method_nr
     694              :       END IF
     695              : 
     696            8 :       SELECT CASE (actual_extrapolation_method_nr)
     697              :       CASE (wfi_linear_wf_method_nr)
     698            8 :          IF (wf_history%snapshot_count < 2) THEN
     699            4 :             actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
     700              :          END IF
     701              :       CASE (wfi_linear_p_method_nr)
     702           12 :          IF (wf_history%snapshot_count < 2) THEN
     703            6 :             actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
     704              :          END IF
     705              :       CASE (wfi_linear_ps_method_nr)
     706        16742 :          IF (wf_history%snapshot_count < 2) THEN
     707            6 :             actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
     708              :          END IF
     709              :       END SELECT
     710              : 
     711        26933 :       IF (PRESENT(extrapolation_method_nr)) THEN
     712        26933 :          extrapolation_method_nr = actual_extrapolation_method_nr
     713              :       END IF
     714        26933 :       my_orthogonal_wf = .FALSE.
     715              : 
     716            8 :       SELECT CASE (actual_extrapolation_method_nr)
     717              :       CASE (wfi_frozen_method_nr)
     718            8 :          CPASSERT(.NOT. do_kpoints)
     719            8 :          t0_state => wfi_get_snapshot(wf_history, wf_index=1)
     720            8 :          CPASSERT(ASSOCIATED(t0_state%rho_frozen))
     721              : 
     722            8 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     723            8 :          CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
     724              : 
     725            8 :          CALL qs_rho_get(t0_state%rho_frozen, rho_ao=rho_frozen_ao)
     726            8 :          CALL qs_rho_get(rho, rho_ao=rho_ao)
     727           16 :          DO ispin = 1, SIZE(rho_frozen_ao)
     728              :             CALL dbcsr_copy(rho_ao(ispin)%matrix, &
     729              :                             rho_frozen_ao(ispin)%matrix, &
     730           16 :                             keep_sparsity=.TRUE.)
     731              :          END DO
     732              :          !FM updating rho_ao directly with t0_state%rho_ao would have the
     733              :          !FM wrong matrix structure
     734            8 :          CALL qs_rho_update_rho(rho, qs_env=qs_env)
     735            8 :          CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
     736              : 
     737            8 :          my_orthogonal_wf = .FALSE.
     738              :       CASE (wfi_use_prev_p_method_nr)
     739          476 :          t0_state => wfi_get_snapshot(wf_history, wf_index=1)
     740          476 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     741          476 :          CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
     742          476 :          IF (do_kpoints) THEN
     743          218 :             CPASSERT(ASSOCIATED(t0_state%rho_ao_kp))
     744          218 :             CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
     745          524 :             DO ispin = 1, SIZE(t0_state%rho_ao_kp, 1)
     746        31248 :                DO img = 1, SIZE(t0_state%rho_ao_kp, 2)
     747        31030 :                   IF (img > SIZE(rho_ao_kp, 2)) THEN
     748              :                      CALL cp_warn(__LOCATION__, &
     749           18 :                                   "Change in cell neighborlist: might affect quality of initial guess")
     750              :                   ELSE
     751              :                      CALL dbcsr_copy(rho_ao_kp(ispin, img)%matrix, &
     752              :                                      t0_state%rho_ao_kp(ispin, img)%matrix, &
     753        30706 :                                      keep_sparsity=.TRUE.)
     754              :                   END IF
     755              :                END DO
     756              :             END DO
     757              :          ELSE
     758          258 :             CPASSERT(ASSOCIATED(t0_state%rho_ao))
     759          258 :             CALL qs_rho_get(rho, rho_ao=rho_ao)
     760          646 :             DO ispin = 1, SIZE(t0_state%rho_ao)
     761              :                CALL dbcsr_copy(rho_ao(ispin)%matrix, &
     762              :                                t0_state%rho_ao(ispin)%matrix, &
     763          646 :                                keep_sparsity=.TRUE.)
     764              :             END DO
     765              :          END IF
     766              :          !FM updating rho_ao directly with t0_state%rho_ao would have the
     767              :          !FM wrong matrix structure
     768          476 :          CALL qs_rho_update_rho(rho, qs_env=qs_env)
     769          476 :          CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
     770              :       CASE (wfi_use_prev_wf_method_nr)
     771          666 :          my_orthogonal_wf = .TRUE.
     772          666 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     773          666 :          CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
     774              : 
     775          666 :          IF (do_kpoints) THEN
     776            6 :             CALL wfi_use_prev_wf_kp(qs_env, io_unit, print_level)
     777              :          ELSE
     778          660 :             CALL qs_rho_get(rho, rho_ao=rho_ao)
     779         1564 :             DO ispin = 1, SIZE(mos)
     780          904 :                CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
     781          904 :                CALL reorthogonalize_vectors(qs_env, v_matrix=mo_coeff, n_col=nmo)
     782         2468 :                CALL calculate_density_matrix(mo_set=mos(ispin), density_matrix=rho_ao(ispin)%matrix)
     783              :             END DO
     784          660 :             CALL qs_rho_update_rho(rho, qs_env=qs_env)
     785          660 :             CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
     786              :          END IF
     787              : 
     788              :       CASE (wfi_use_guess_method_nr)
     789              :          !FM more clean to do it here, but it
     790              :          !FM might need to read a file (restart) and thus globenv
     791              :          !FM I do not want globenv here, thus done by the caller
     792              :          !FM (btw. it also needs the eigensolver, and unless you relocate it
     793              :          !FM gives circular dependencies)
     794        10337 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     795        10337 :          CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
     796              :       CASE (wfi_linear_wf_method_nr)
     797            4 :          CPASSERT(.NOT. do_kpoints)
     798            4 :          t0_state => wfi_get_snapshot(wf_history, wf_index=2)
     799            4 :          t1_state => wfi_get_snapshot(wf_history, wf_index=1)
     800            4 :          CPASSERT(ASSOCIATED(t0_state))
     801            4 :          CPASSERT(ASSOCIATED(t1_state))
     802            4 :          CPASSERT(ASSOCIATED(t0_state%wf))
     803            4 :          CPASSERT(ASSOCIATED(t1_state%wf))
     804            4 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     805            4 :          CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
     806              : 
     807            4 :          my_orthogonal_wf = .TRUE.
     808            4 :          t0 = 0.0_dp
     809            4 :          t1 = t1_state%dt
     810            4 :          t2 = t1 + dt
     811            4 :          CALL qs_rho_get(rho, rho_ao=rho_ao)
     812            8 :          DO ispin = 1, SIZE(mos)
     813              :             CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, &
     814            4 :                             nmo=nmo)
     815              :             CALL cp_fm_scale_and_add(alpha=0.0_dp, &
     816              :                                      matrix_a=mo_coeff, &
     817              :                                      matrix_b=t1_state%wf(ispin), &
     818            4 :                                      beta=(t2 - t0)/(t1 - t0))
     819              :             ! this copy should be unnecessary
     820              :             CALL cp_fm_scale_and_add(alpha=1.0_dp, &
     821              :                                      matrix_a=mo_coeff, &
     822            4 :                                      beta=(t1 - t2)/(t1 - t0), matrix_b=t0_state%wf(ispin))
     823              :             CALL reorthogonalize_vectors(qs_env, &
     824              :                                          v_matrix=mo_coeff, &
     825            4 :                                          n_col=nmo)
     826              :             CALL calculate_density_matrix(mo_set=mos(ispin), &
     827           12 :                                           density_matrix=rho_ao(ispin)%matrix)
     828              :          END DO
     829            4 :          CALL qs_rho_update_rho(rho, qs_env=qs_env)
     830              : 
     831              :          CALL qs_ks_did_change(qs_env%ks_env, &
     832            4 :                                rho_changed=.TRUE.)
     833              :       CASE (wfi_linear_p_method_nr)
     834            6 :          t0_state => wfi_get_snapshot(wf_history, wf_index=2)
     835            6 :          t1_state => wfi_get_snapshot(wf_history, wf_index=1)
     836            6 :          CPASSERT(ASSOCIATED(t0_state))
     837            6 :          CPASSERT(ASSOCIATED(t1_state))
     838            6 :          IF (do_kpoints) THEN
     839            2 :             CPASSERT(ASSOCIATED(t0_state%rho_ao_kp))
     840            2 :             CPASSERT(ASSOCIATED(t1_state%rho_ao_kp))
     841              :          ELSE
     842            4 :             CPASSERT(ASSOCIATED(t0_state%rho_ao))
     843            4 :             CPASSERT(ASSOCIATED(t1_state%rho_ao))
     844              :          END IF
     845            6 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     846            6 :          CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
     847              : 
     848            6 :          t0 = 0.0_dp
     849            6 :          t1 = t1_state%dt
     850            6 :          t2 = t1 + dt
     851            6 :          IF (do_kpoints) THEN
     852            2 :             CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
     853            4 :             DO ispin = 1, SIZE(rho_ao_kp, 1)
     854          528 :                DO img = 1, SIZE(rho_ao_kp, 2)
     855          524 :                   IF (img > SIZE(t0_state%rho_ao_kp, 2) .OR. &
     856            2 :                       img > SIZE(t1_state%rho_ao_kp, 2)) THEN
     857              :                      CALL cp_warn(__LOCATION__, &
     858           22 :                                   "Change in cell neighborlist: might affect quality of initial guess")
     859              :                   ELSE
     860              :                      CALL dbcsr_add(rho_ao_kp(ispin, img)%matrix, t1_state%rho_ao_kp(ispin, img)%matrix, &
     861          502 :                                     alpha_scalar=0.0_dp, beta_scalar=(t2 - t0)/(t1 - t0)) ! this copy should be unnecessary
     862              :                      CALL dbcsr_add(rho_ao_kp(ispin, img)%matrix, t0_state%rho_ao_kp(ispin, img)%matrix, &
     863          502 :                                     alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0))
     864              :                   END IF
     865              :                END DO
     866              :             END DO
     867              :          ELSE
     868            4 :             CALL qs_rho_get(rho, rho_ao=rho_ao)
     869            8 :             DO ispin = 1, SIZE(rho_ao)
     870              :                CALL dbcsr_add(rho_ao(ispin)%matrix, t1_state%rho_ao(ispin)%matrix, &
     871            4 :                               alpha_scalar=0.0_dp, beta_scalar=(t2 - t0)/(t1 - t0)) ! this copy should be unnecessary
     872              :                CALL dbcsr_add(rho_ao(ispin)%matrix, t0_state%rho_ao(ispin)%matrix, &
     873            8 :                               alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0))
     874              :             END DO
     875              :          END IF
     876            6 :          CALL qs_rho_update_rho(rho, qs_env=qs_env)
     877            6 :          CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
     878              :       CASE (wfi_linear_ps_method_nr)
     879              :          ! wf not calculated, extract with PSC renormalized?
     880              :          ! use wf_linear?
     881           12 :          CPASSERT(.NOT. do_kpoints)
     882           12 :          t0_state => wfi_get_snapshot(wf_history, wf_index=2)
     883           12 :          t1_state => wfi_get_snapshot(wf_history, wf_index=1)
     884           12 :          CPASSERT(ASSOCIATED(t0_state))
     885           12 :          CPASSERT(ASSOCIATED(t1_state))
     886           12 :          CPASSERT(ASSOCIATED(t0_state%wf))
     887           12 :          CPASSERT(ASSOCIATED(t1_state%wf))
     888           12 :          IF (wf_history%store_overlap) THEN
     889            4 :             CPASSERT(ASSOCIATED(t0_state%overlap))
     890            4 :             CPASSERT(ASSOCIATED(t1_state%overlap))
     891              :          END IF
     892           12 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     893           12 :          IF (nvec >= wf_history%memory_depth) THEN
     894           12 :             IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. (qs_env%scf_control%eps_scf_hist /= 0)) THEN
     895            0 :                qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
     896            0 :                qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
     897            0 :                qs_env%scf_control%outer_scf%have_scf = .FALSE.
     898           12 :             ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
     899            0 :                qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
     900            0 :                qs_env%scf_control%outer_scf%have_scf = .FALSE.
     901           12 :             ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
     902            0 :                qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
     903              :             END IF
     904              :          END IF
     905              : 
     906           12 :          my_orthogonal_wf = .TRUE.
     907              :          ! use PS_2=2 PS_1-PS_0
     908              :          ! C_2 comes from using PS_2 as a projector acting on C_1
     909           12 :          CALL qs_rho_get(rho, rho_ao=rho_ao)
     910           24 :          DO ispin = 1, SIZE(mos)
     911           12 :             NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
     912           12 :             CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
     913              :             CALL cp_fm_get_info(mo_coeff, nrow_global=n, ncol_global=k, &
     914           12 :                                 matrix_struct=matrix_struct)
     915              :             CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, &
     916           12 :                                      nrow_global=k, ncol_global=k)
     917           12 :             CALL cp_fm_create(csc, matrix_struct_new)
     918           12 :             CALL cp_fm_struct_release(matrix_struct_new)
     919              : 
     920           12 :             IF (use_overlap) THEN
     921            4 :                CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), mo_coeff, k)
     922            4 :                CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), mo_coeff, 0.0_dp, csc)
     923              :             ELSE
     924              :                CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
     925            8 :                                   t1_state%wf(ispin), 0.0_dp, csc)
     926              :             END IF
     927           12 :             CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, mo_coeff)
     928           12 :             CALL cp_fm_release(csc)
     929           12 :             CALL cp_fm_scale_and_add(-1.0_dp, mo_coeff, 2.0_dp, t1_state%wf(ispin))
     930              :             CALL reorthogonalize_vectors(qs_env, &
     931              :                                          v_matrix=mo_coeff, &
     932           12 :                                          n_col=k)
     933              :             CALL calculate_density_matrix(mo_set=mos(ispin), &
     934           48 :                                           density_matrix=rho_ao(ispin)%matrix)
     935              :          END DO
     936           12 :          CALL qs_rho_update_rho(rho, qs_env=qs_env)
     937           12 :          CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
     938              : 
     939              :       CASE (wfi_ps_method_nr)
     940              :          ! figure out the actual number of vectors to use in the extrapolation:
     941          376 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     942          376 :          CPASSERT(nvec > 0)
     943          376 :          IF (nvec >= wf_history%memory_depth) THEN
     944          178 :             IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. (qs_env%scf_control%eps_scf_hist /= 0)) THEN
     945            0 :                qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
     946            0 :                qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
     947            0 :                qs_env%scf_control%outer_scf%have_scf = .FALSE.
     948          178 :             ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
     949            0 :                qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
     950            0 :                qs_env%scf_control%outer_scf%have_scf = .FALSE.
     951          178 :             ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
     952            0 :                qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
     953              :             END IF
     954              :          END IF
     955              : 
     956          376 :          IF (do_kpoints) THEN
     957            4 :             CALL wfi_extrapolate_ps_aspc_kp(wf_history, qs_env, nvec, io_unit, print_level)
     958            4 :             my_orthogonal_wf = .TRUE.
     959              :          ELSE
     960          372 :             my_orthogonal_wf = .TRUE.
     961          822 :             DO ispin = 1, SIZE(mos)
     962          450 :                NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
     963          450 :                CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
     964              :                CALL cp_fm_get_info(mo_coeff, nrow_global=n, ncol_global=k, &
     965          450 :                                    matrix_struct=matrix_struct)
     966          450 :                CALL cp_fm_create(fm_tmp, matrix_struct)
     967              :                CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, &
     968          450 :                                         nrow_global=k, ncol_global=k)
     969          450 :                CALL cp_fm_create(csc, matrix_struct_new)
     970          450 :                CALL cp_fm_struct_release(matrix_struct_new)
     971              :                ! first the most recent
     972          450 :                t1_state => wfi_get_snapshot(wf_history, wf_index=1)
     973          450 :                CALL cp_fm_to_fm(t1_state%wf(ispin), mo_coeff)
     974          450 :                alpha = nvec
     975          450 :                CALL cp_fm_scale(alpha, mo_coeff)
     976          450 :                CALL qs_rho_get(rho, rho_ao=rho_ao)
     977          962 :                DO i = 2, nvec
     978          512 :                   t0_state => wfi_get_snapshot(wf_history, wf_index=i)
     979          512 :                   IF (use_overlap) THEN
     980          474 :                      CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
     981          474 :                      CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
     982              :                   ELSE
     983              :                      CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
     984           38 :                                         t1_state%wf(ispin), 0.0_dp, csc)
     985              :                   END IF
     986          512 :                   CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
     987          512 :                   alpha = -1.0_dp*alpha*REAL(nvec - i + 1, dp)/REAL(i, dp)
     988          962 :                   CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp)
     989              :                END DO
     990              : 
     991          450 :                CALL cp_fm_release(csc)
     992          450 :                CALL cp_fm_release(fm_tmp)
     993              :                CALL reorthogonalize_vectors(qs_env, &
     994              :                                             v_matrix=mo_coeff, &
     995          450 :                                             n_col=k)
     996              :                CALL calculate_density_matrix(mo_set=mos(ispin), &
     997         1722 :                                              density_matrix=rho_ao(ispin)%matrix)
     998              :             END DO
     999          372 :             CALL qs_rho_update_rho(rho, qs_env=qs_env)
    1000          372 :             CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
    1001              :          END IF
    1002              : 
    1003              :       CASE (wfi_aspc_nr)
    1004        14874 :          CALL cite_reference(Kolafa2004)
    1005        14874 :          CALL cite_reference(Kuhne2007)
    1006              :          ! figure out the actual number of vectors to use in the extrapolation:
    1007        14874 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
    1008        14874 :          CPASSERT(nvec > 0)
    1009        14874 :          IF (nvec >= wf_history%memory_depth) THEN
    1010         9652 :             IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. &
    1011              :                 (qs_env%scf_control%eps_scf_hist /= 0)) THEN
    1012           18 :                qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
    1013           18 :                qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
    1014           18 :                qs_env%scf_control%outer_scf%have_scf = .FALSE.
    1015         9634 :             ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
    1016           62 :                qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
    1017           62 :                qs_env%scf_control%outer_scf%have_scf = .FALSE.
    1018         9572 :             ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
    1019            8 :                qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
    1020              :             END IF
    1021              :          END IF
    1022              : 
    1023        14874 :          IF (do_kpoints) THEN
    1024          420 :             CALL wfi_extrapolate_ps_aspc_kp(wf_history, qs_env, nvec, io_unit, print_level)
    1025          420 :             my_orthogonal_wf = .TRUE.
    1026              :          ELSE
    1027        14454 :             my_orthogonal_wf = .TRUE.
    1028        14454 :             CALL qs_rho_get(rho, rho_ao=rho_ao)
    1029        30127 :             DO ispin = 1, SIZE(mos)
    1030        15673 :                NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
    1031        15673 :                CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
    1032              :                CALL cp_fm_get_info(mo_coeff, &
    1033              :                                    nrow_global=n, &
    1034              :                                    ncol_global=k, &
    1035        15673 :                                    matrix_struct=matrix_struct)
    1036        15673 :                CALL cp_fm_create(fm_tmp, matrix_struct, set_zero=.TRUE.)
    1037              :                CALL cp_fm_struct_create(matrix_struct_new, &
    1038              :                                         template_fmstruct=matrix_struct, &
    1039              :                                         nrow_global=k, &
    1040        15673 :                                         ncol_global=k)
    1041        15673 :                CALL cp_fm_create(csc, matrix_struct_new, set_zero=.TRUE.)
    1042        15673 :                CALL cp_fm_struct_release(matrix_struct_new)
    1043              :                ! first the most recent
    1044              :                t1_state => wfi_get_snapshot(wf_history, &
    1045        15673 :                                             wf_index=1)
    1046        15673 :                CALL cp_fm_to_fm(t1_state%wf(ispin), mo_coeff)
    1047        15673 :                alpha = REAL(4*nvec - 2, KIND=dp)/REAL(nvec + 1, KIND=dp)
    1048        15673 :                IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
    1049              :                   WRITE (UNIT=io_unit, FMT="(/,T2,A,/,/,T3,A,I0,/,/,T3,A2,I0,A4,F10.6)") &
    1050         3138 :                      "Parameters for the always stable predictor-corrector (ASPC) method:", &
    1051         3138 :                      "ASPC order: ", MAX(nvec - 2, 0), &
    1052         6276 :                      "B(", 1, ") = ", alpha
    1053              :                END IF
    1054        15673 :                CALL cp_fm_scale(alpha, mo_coeff)
    1055              : 
    1056        61229 :                DO i = 2, nvec
    1057        45556 :                   t0_state => wfi_get_snapshot(wf_history, wf_index=i)
    1058        45556 :                   IF (use_overlap) THEN
    1059        34016 :                      CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
    1060        34016 :                      CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
    1061              :                   ELSE
    1062              :                      CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
    1063        11540 :                                         t1_state%wf(ispin), 0.0_dp, csc)
    1064              :                   END IF
    1065        45556 :                   CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
    1066              :                   alpha = (-1.0_dp)**(i + 1)*REAL(i, KIND=dp)* &
    1067        45556 :                           binomial(2*nvec, nvec - i)/binomial(2*nvec - 2, nvec - 1)
    1068        45556 :                   IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
    1069              :                      WRITE (UNIT=io_unit, FMT="(T3,A2,I0,A4,F10.6)") &
    1070         9414 :                         "B(", i, ") = ", alpha
    1071              :                   END IF
    1072        61229 :                   CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp)
    1073              :                END DO
    1074        15673 :                CALL cp_fm_release(csc)
    1075        15673 :                CALL cp_fm_release(fm_tmp)
    1076              :                CALL reorthogonalize_vectors(qs_env, &
    1077              :                                             v_matrix=mo_coeff, &
    1078        15673 :                                             n_col=k)
    1079              :                CALL calculate_density_matrix(mo_set=mos(ispin), &
    1080        45800 :                                              density_matrix=rho_ao(ispin)%matrix)
    1081              :             END DO
    1082        14454 :             CALL qs_rho_update_rho(rho, qs_env=qs_env)
    1083        14454 :             CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
    1084              :          END IF ! do_kpoints
    1085              : 
    1086              :       CASE (wfi_gext_proj_nr)
    1087          146 :          IF (do_kpoints) THEN
    1088           10 :             nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
    1089           10 :             CPASSERT(nvec > 0)
    1090           10 :             CALL wfi_extrapolate_gext_proj_kp(wf_history, qs_env, nvec, io_unit, print_level)
    1091           10 :             my_orthogonal_wf = .TRUE.
    1092              :          ELSE
    1093              : 
    1094              :             ! figure out the actual number of vectors to use in the extrapolation:
    1095          136 :             nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
    1096          136 :             IF (nvec >= wf_history%memory_depth) THEN
    1097           88 :                IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. &
    1098              :                    (qs_env%scf_control%eps_scf_hist /= 0)) THEN
    1099            0 :                   qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
    1100            0 :                   qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
    1101            0 :                   qs_env%scf_control%outer_scf%have_scf = .FALSE.
    1102           88 :                ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
    1103            0 :                   qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
    1104            0 :                   qs_env%scf_control%outer_scf%have_scf = .FALSE.
    1105           88 :                ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
    1106            0 :                   qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
    1107              :                END IF
    1108              :             END IF
    1109          136 :             CPASSERT(nvec > 0)
    1110              : 
    1111              :             ! get the coefficients for the fitting
    1112          408 :             ALLOCATE (coeffs(nvec))
    1113          136 :             NULLIFY (matrix_s)
    1114          136 :             CALL get_qs_env(qs_env, matrix_s=matrix_s)
    1115              :             CALL diff_fitting(wf_history, matrix_s(1)%matrix, coeffs, nvec, &
    1116          136 :                               1e-4_dp, io_unit, print_level)
    1117              : 
    1118          136 :             my_orthogonal_wf = .TRUE.
    1119          136 :             CALL qs_rho_get(rho, rho_ao=rho_ao)
    1120          328 :             DO ispin = 1, SIZE(mos)
    1121          192 :                NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
    1122          192 :                CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
    1123              :                CALL cp_fm_get_info(mo_coeff, &
    1124              :                                    nrow_global=n, &
    1125              :                                    ncol_global=k, &
    1126          192 :                                    matrix_struct=matrix_struct)
    1127          192 :                CALL cp_fm_create(fm_tmp, matrix_struct)
    1128              :                CALL cp_fm_struct_create(matrix_struct_new, &
    1129              :                                         template_fmstruct=matrix_struct, &
    1130              :                                         nrow_global=k, &
    1131          192 :                                         ncol_global=k)
    1132          192 :                CALL cp_fm_create(csc, matrix_struct_new)
    1133          192 :                CALL cp_fm_struct_release(matrix_struct_new)
    1134              : 
    1135          192 :                t1_state => wfi_get_snapshot(wf_history, wf_index=1)
    1136              : 
    1137              :                ! do the linear combination of previous PSs
    1138          192 :                CALL cp_fm_set_all(mo_coeff, 0.0_dp)
    1139          704 :                DO i = 1, nvec
    1140          512 :                   t0_state => wfi_get_snapshot(wf_history, wf_index=i)
    1141          512 :                   CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
    1142          512 :                   CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
    1143          512 :                   CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
    1144          704 :                   CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, coeffs(i), fm_tmp)
    1145              :                END DO
    1146          192 :                CALL cp_fm_release(csc)
    1147          192 :                CALL cp_fm_release(fm_tmp)
    1148              :                CALL reorthogonalize_vectors(qs_env, &
    1149              :                                             v_matrix=mo_coeff, &
    1150          192 :                                             n_col=k)
    1151              :                CALL calculate_density_matrix(mo_set=mos(ispin), &
    1152          712 :                                              density_matrix=rho_ao(ispin)%matrix)
    1153              :             END DO
    1154          136 :             CALL qs_rho_update_rho(rho, qs_env=qs_env)
    1155          136 :             CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
    1156              : 
    1157          136 :             DEALLOCATE (coeffs)
    1158              : 
    1159              :          END IF
    1160              : 
    1161              :       CASE (wfi_gext_proj_qtr_nr)
    1162           28 :          IF (do_kpoints) THEN
    1163            4 :             nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
    1164            4 :             CPASSERT(nvec > 0)
    1165            4 :             CALL wfi_extrapolate_gext_proj_kp(wf_history, qs_env, nvec, io_unit, print_level)
    1166            4 :             my_orthogonal_wf = .TRUE.
    1167              :          ELSE
    1168              : 
    1169              :             ! figure out the actual number of vectors to use in the extrapolation:
    1170           24 :             nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
    1171           24 :             IF (nvec >= wf_history%memory_depth) THEN
    1172            8 :                IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. &
    1173              :                    (qs_env%scf_control%eps_scf_hist /= 0)) THEN
    1174            0 :                   qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
    1175            0 :                   qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
    1176            0 :                   qs_env%scf_control%outer_scf%have_scf = .FALSE.
    1177            8 :                ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
    1178            0 :                   qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
    1179            0 :                   qs_env%scf_control%outer_scf%have_scf = .FALSE.
    1180            8 :                ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
    1181            0 :                   qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
    1182              :                END IF
    1183              :             END IF
    1184           24 :             CPASSERT(nvec > 0)
    1185              : 
    1186              :             ! get the coefficients for the fitting
    1187           72 :             ALLOCATE (coeffs(nvec))
    1188           24 :             NULLIFY (matrix_s)
    1189           24 :             CALL get_qs_env(qs_env, matrix_s=matrix_s)
    1190              :             CALL tr_fitting(wf_history, matrix_s(1)%matrix, coeffs, nvec, &
    1191           24 :                             1e-4_dp, io_unit, print_level)
    1192              : 
    1193           24 :             my_orthogonal_wf = .TRUE.
    1194           24 :             CALL qs_rho_get(rho, rho_ao=rho_ao)
    1195           48 :             DO ispin = 1, SIZE(mos)
    1196           24 :                NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
    1197           24 :                CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
    1198              :                CALL cp_fm_get_info(mo_coeff, &
    1199              :                                    nrow_global=n, &
    1200              :                                    ncol_global=k, &
    1201           24 :                                    matrix_struct=matrix_struct)
    1202           24 :                CALL cp_fm_create(fm_tmp, matrix_struct)
    1203              :                CALL cp_fm_struct_create(matrix_struct_new, &
    1204              :                                         template_fmstruct=matrix_struct, &
    1205              :                                         nrow_global=k, &
    1206           24 :                                         ncol_global=k)
    1207           24 :                CALL cp_fm_create(csc, matrix_struct_new)
    1208           24 :                CALL cp_fm_struct_release(matrix_struct_new)
    1209              : 
    1210           24 :                t1_state => wfi_get_snapshot(wf_history, wf_index=1)
    1211              : 
    1212              :                ! do the linear combination of previous PSs
    1213           24 :                CALL cp_fm_set_all(mo_coeff, 0.0_dp)
    1214          104 :                DO i = 1, nvec
    1215           80 :                   t0_state => wfi_get_snapshot(wf_history, wf_index=i)
    1216           80 :                   CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
    1217           80 :                   CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
    1218           80 :                   CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
    1219          104 :                   CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, coeffs(i), fm_tmp)
    1220              :                END DO
    1221           24 :                CALL cp_fm_release(csc)
    1222           24 :                CALL cp_fm_release(fm_tmp)
    1223              :                CALL reorthogonalize_vectors(qs_env, &
    1224              :                                             v_matrix=mo_coeff, &
    1225           24 :                                             n_col=k)
    1226              :                CALL calculate_density_matrix(mo_set=mos(ispin), &
    1227           96 :                                              density_matrix=rho_ao(ispin)%matrix)
    1228              :             END DO
    1229           24 :             CALL qs_rho_update_rho(rho, qs_env=qs_env)
    1230           24 :             CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
    1231              : 
    1232           24 :             DEALLOCATE (coeffs)
    1233              : 
    1234              :          END IF
    1235              : 
    1236              :       CASE default
    1237              :          CALL cp_abort(__LOCATION__, &
    1238              :                        "Unknown interpolation method: "// &
    1239        26933 :                        TRIM(ADJUSTL(cp_to_string(wf_history%interpolation_method_nr))))
    1240              :       END SELECT
    1241        26933 :       IF (PRESENT(orthogonal_wf)) orthogonal_wf = my_orthogonal_wf
    1242              :       CALL cp_print_key_finished_output(io_unit, logger, qs_env%input, &
    1243        26933 :                                         "DFT%SCF%PRINT%PROGRAM_RUN_INFO")
    1244        26933 :       CALL timestop(handle)
    1245        26933 :    END SUBROUTINE wfi_extrapolate
    1246              : 
    1247              : ! **************************************************************************************************
    1248              : !> \brief Reorthogonalizes the wavefunctions from the previous step for k-points
    1249              : !>        using the current S(k) metric and rebuilds the density matrix.
    1250              : !> \param qs_env The QS environment
    1251              : !> \param io_unit output unit
    1252              : !> \param print_level print level
    1253              : !> \param pbc_shift_ref ...
    1254              : !> \param load_snapshot_wf ...
    1255              : ! **************************************************************************************************
    1256          444 :    SUBROUTINE wfi_use_prev_wf_kp(qs_env, io_unit, print_level, pbc_shift_ref, load_snapshot_wf)
    1257              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1258              :       INTEGER, INTENT(IN)                                :: io_unit, print_level
    1259              :       INTEGER, DIMENSION(:, :), INTENT(IN), OPTIONAL     :: pbc_shift_ref
    1260              :       LOGICAL, INTENT(IN), OPTIONAL                      :: load_snapshot_wf
    1261              : 
    1262              :       CHARACTER(len=*), PARAMETER :: routineN = 'wfi_use_prev_wf_kp'
    1263              : 
    1264          444 :       COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:)        :: col_scaling
    1265              :       INTEGER                                            :: chol_info, handle, ik, ikp, ispin, j, &
    1266              :                                                             kplocal, nao, nmo, nspin
    1267          444 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: pbc_shift_cur, pbc_shift_src
    1268              :       INTEGER, DIMENSION(2)                              :: kp_range
    1269              :       LOGICAL                                            :: reload_snapshot_wf, use_pbc_phase_ref, &
    1270              :                                                             use_real_wfn
    1271              :       REAL(KIND=dp)                                      :: eval_thresh
    1272          444 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eigenvalues
    1273          444 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: xkp
    1274              :       TYPE(cp_cfm_type)                                  :: cfm_evecs, cfm_mhalf, cfm_nao_nmo_work, &
    1275              :                                                             cmos_new, csc_cfm
    1276          444 :       TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:)       :: csmat_cur
    1277          444 :       TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER     :: ao_ao_fm_pools_kp
    1278              :       TYPE(cp_fm_struct_type), POINTER                   :: ao_ao_struct, nmo_nmo_struct
    1279              :       TYPE(cp_fm_type), POINTER                          :: imos, mo_coeff, rmos
    1280          444 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s_kp
    1281              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1282              :       TYPE(kpoint_env_type), POINTER                     :: kp
    1283              :       TYPE(kpoint_type), POINTER                         :: kpoints
    1284              :       TYPE(qs_matrix_pools_type), POINTER                :: mpools_kp
    1285              :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
    1286              :       TYPE(qs_wf_history_type), POINTER                  :: wf_history
    1287              :       TYPE(qs_wf_snapshot_type), POINTER                 :: t1_state
    1288              :       TYPE(scf_control_type), POINTER                    :: scf_control
    1289              : 
    1290          444 :       CALL timeset(routineN, handle)
    1291              : 
    1292          444 :       NULLIFY (dft_control, scf_control, kpoints, matrix_s_kp, scf_env, kp, &
    1293          444 :                mo_coeff, rmos, imos, wf_history, t1_state)
    1294              : 
    1295              :       CALL get_qs_env(qs_env, kpoints=kpoints, matrix_s_kp=matrix_s_kp, scf_env=scf_env, &
    1296          444 :                       dft_control=dft_control, scf_control=scf_control)
    1297          444 :       CALL get_kpoint_info(kpoints, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range)
    1298          444 :       kplocal = kp_range(2) - kp_range(1) + 1
    1299              : 
    1300          444 :       IF (use_real_wfn) THEN
    1301            0 :          CALL timestop(handle)
    1302            0 :          RETURN
    1303              :       END IF
    1304              : 
    1305          444 :       wf_history => qs_env%wf_history
    1306          444 :       reload_snapshot_wf = .FALSE.
    1307          444 :       IF (PRESENT(load_snapshot_wf)) reload_snapshot_wf = load_snapshot_wf
    1308          444 :       IF (PRESENT(pbc_shift_ref)) THEN
    1309         1314 :          ALLOCATE (pbc_shift_src(3, SIZE(pbc_shift_ref, 2)))
    1310        12086 :          pbc_shift_src(:, :) = pbc_shift_ref(:, :)
    1311          442 :          use_pbc_phase_ref = .TRUE.
    1312              :       ELSE
    1313            6 :          use_pbc_phase_ref = .FALSE.
    1314            6 :          IF (ASSOCIATED(wf_history)) THEN
    1315            6 :             IF (wf_history%store_wf_kp .AND. wf_history%snapshot_count > 0) THEN
    1316            4 :                t1_state => wfi_get_snapshot(wf_history, wf_index=1)
    1317            4 :                CPASSERT(ASSOCIATED(t1_state%wf_kp))
    1318            4 :                CPASSERT(ASSOCIATED(t1_state%kp_pbc_shift))
    1319            4 :                reload_snapshot_wf = .TRUE.
    1320           12 :                ALLOCATE (pbc_shift_src(3, SIZE(t1_state%kp_pbc_shift, 2)))
    1321          132 :                pbc_shift_src(:, :) = t1_state%kp_pbc_shift(:, :)
    1322              :                use_pbc_phase_ref = .TRUE.
    1323              :             END IF
    1324              :          END IF
    1325              :       END IF
    1326          442 :       IF (use_pbc_phase_ref) CALL wfi_compute_kp_pbc_shift(qs_env, pbc_shift_cur)
    1327              : 
    1328          444 :       kp => kpoints%kp_env(1)%kpoint_env
    1329          444 :       nspin = SIZE(kp%mos, 2)
    1330          444 :       CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo, mo_coeff=mo_coeff)
    1331              : 
    1332          444 :       IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
    1333              :          WRITE (UNIT=io_unit, FMT="(/,T2,A)") &
    1334            0 :             "Using previous wavefunctions as initial guess for k-points (with reorthogonalization)"
    1335              :       END IF
    1336              : 
    1337          444 :       CALL get_kpoint_info(kpoints, mpools=mpools_kp)
    1338          444 :       CALL mpools_get(mpools_kp, ao_ao_fm_pools=ao_ao_fm_pools_kp)
    1339          444 :       ao_ao_struct => fm_pool_get_el_struct(ao_ao_fm_pools_kp(1)%pool)
    1340              : 
    1341          444 :       CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct)
    1342          444 :       CALL cp_cfm_create(cfm_nao_nmo_work, mo_coeff%matrix_struct)
    1343              : 
    1344          444 :       NULLIFY (nmo_nmo_struct)
    1345              :       CALL cp_fm_struct_create(nmo_nmo_struct, template_fmstruct=mo_coeff%matrix_struct, &
    1346          444 :                                nrow_global=nmo, ncol_global=nmo)
    1347          444 :       CALL cp_cfm_create(csc_cfm, nmo_nmo_struct)
    1348          444 :       CALL cp_fm_struct_release(nmo_nmo_struct)
    1349              : 
    1350         2420 :       ALLOCATE (csmat_cur(kplocal))
    1351         1532 :       DO ikp = 1, kplocal
    1352         1532 :          CALL cp_cfm_create(csmat_cur(ikp), ao_ao_struct)
    1353              :       END DO
    1354              : 
    1355              :       CALL kpoint_operator_get_all(kpoints, scf_env%scf_work1(1), matrix_s_kp, kpoint_slot_s, &
    1356          444 :                                    csmat_cur, matrix_row=kpoint_spin_free)
    1357              : 
    1358              :       ! Phase B: bring the WFN from its saved/internal PBC image convention to
    1359              :       ! the current convention, then orthogonalize it with respect to S(k).
    1360         1332 :       ALLOCATE (eigenvalues(nmo))
    1361          444 :       eval_thresh = 1.0E-12_dp
    1362              : 
    1363         1532 :       DO ikp = 1, kplocal
    1364         1088 :          kp => kpoints%kp_env(ikp)%kpoint_env
    1365         2376 :          DO ispin = 1, nspin
    1366         1288 :             CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
    1367         1288 :             CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
    1368         1288 :             IF (reload_snapshot_wf) THEN
    1369           16 :                CALL cp_fm_to_fm(t1_state%wf_kp(ikp, 1, ispin), rmos)
    1370           16 :                CALL cp_fm_to_fm(t1_state%wf_kp(ikp, 2, ispin), imos)
    1371              :             END IF
    1372         1288 :             IF (use_pbc_phase_ref) THEN
    1373         1280 :                ik = kp_range(1) + ikp - 1
    1374              :                CALL wfi_apply_kp_pbc_phase_fm(rmos, imos, pbc_shift_cur - pbc_shift_src, &
    1375        25224 :                                               xkp(1:3, ik), matrix_s_kp(1, 1)%matrix)
    1376              :             END IF
    1377         1288 :             CALL cp_fm_to_cfm(rmos, imos, cmos_new)
    1378              : 
    1379              :             CALL cp_cfm_gemm('N', 'N', nao, nmo, nao, z_one, &
    1380         1288 :                              csmat_cur(ikp), cmos_new, z_zero, cfm_nao_nmo_work)
    1381              :             CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, &
    1382         1288 :                              cmos_new, cfm_nao_nmo_work, z_zero, csc_cfm)
    1383              : 
    1384         1288 :             CALL cp_cfm_cholesky_decompose(csc_cfm, info_out=chol_info)
    1385         1288 :             IF (chol_info == 0) THEN
    1386         1284 :                CALL cp_cfm_triangular_multiply(csc_cfm, cmos_new, side='R', invert_tr=.TRUE., uplo_tr='U')
    1387              :             ELSE
    1388            4 :                CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, cmos_new, cfm_nao_nmo_work, z_zero, csc_cfm)
    1389            4 :                CALL cp_cfm_create(cfm_evecs, csc_cfm%matrix_struct)
    1390            4 :                CALL cp_cfm_create(cfm_mhalf, csc_cfm%matrix_struct)
    1391            4 :                CALL cp_cfm_heevd(csc_cfm, cfm_evecs, eigenvalues)
    1392            4 :                CALL cp_cfm_to_cfm(cfm_evecs, cfm_mhalf)
    1393           12 :                ALLOCATE (col_scaling(nmo))
    1394           12 :                DO j = 1, nmo
    1395           12 :                   IF (eigenvalues(j) > eval_thresh) THEN
    1396            4 :                      col_scaling(j) = CMPLX(1.0_dp/SQRT(eigenvalues(j)), 0.0_dp, KIND=dp)
    1397              :                   ELSE
    1398            4 :                      col_scaling(j) = z_zero
    1399              :                   END IF
    1400              :                END DO
    1401            4 :                CALL cp_cfm_column_scale(cfm_mhalf, col_scaling)
    1402            4 :                DEALLOCATE (col_scaling)
    1403            4 :                CALL cp_cfm_gemm('N', 'C', nmo, nmo, nmo, z_one, cfm_mhalf, cfm_evecs, z_zero, csc_cfm)
    1404            4 :                CALL cp_cfm_gemm('N', 'N', nao, nmo, nmo, z_one, cmos_new, csc_cfm, z_zero, cfm_nao_nmo_work)
    1405            4 :                CALL cp_cfm_to_cfm(cfm_nao_nmo_work, cmos_new)
    1406            4 :                CALL cp_cfm_release(cfm_evecs)
    1407            8 :                CALL cp_cfm_release(cfm_mhalf)
    1408              :             END IF
    1409         3664 :             CALL cp_cfm_to_fm(cmos_new, rmos, imos)
    1410              :          END DO
    1411              :          ! the MOS now hold extrapolated coefficients that are orthonormal under
    1412              :          ! the S(k) of the current geometry: a valid trial subspace for solvers.
    1413              :          ! All WFN-based k-point extrapolation methods converge here (PS/ASPC
    1414              :          ! and GEXT_PROJ finish through this routine). The flag is set once
    1415              :          ! at their common write-back.
    1416         1532 :          kp%mos_prefilled = .TRUE.
    1417              :       END DO
    1418          444 :       DEALLOCATE (eigenvalues)
    1419              : 
    1420              :       ! Phase C: Rebuild Density Matrix P(R)
    1421              :       ! Fixed-occupation OT must preserve its occupied rank after extrapolation; its energy
    1422              :       ! labels and any padded UKS columns are not a spectrum for refilling the bands.
    1423              :       CALL qs_kpoint_state_commit(qs_env, update_occupations=.TRUE., &
    1424              :                                   fixed_occupations=scf_env%method == ot_method_nr .AND. &
    1425          884 :                                   .NOT. (dft_control%smear .OR. scf_control%smear%do_smear))
    1426              : 
    1427              :       ! Cleanup
    1428         1532 :       DO ikp = 1, kplocal
    1429         1532 :          CALL cp_cfm_release(csmat_cur(ikp))
    1430              :       END DO
    1431          444 :       DEALLOCATE (csmat_cur)
    1432          444 :       CALL cp_cfm_release(cmos_new)
    1433          444 :       CALL cp_cfm_release(cfm_nao_nmo_work)
    1434          444 :       CALL cp_cfm_release(csc_cfm)
    1435          444 :       IF (ALLOCATED(pbc_shift_cur)) DEALLOCATE (pbc_shift_cur)
    1436          444 :       IF (ALLOCATED(pbc_shift_src)) DEALLOCATE (pbc_shift_src)
    1437              : 
    1438          444 :       CALL timestop(handle)
    1439         2220 :    END SUBROUTINE wfi_use_prev_wf_kp
    1440              : 
    1441              : ! **************************************************************************************************
    1442              : !> \brief Stores the internal PBC image shift used for k-point neighbor-list construction.
    1443              : !>        shift = scaled(pbc(r))-scaled(r), i.e. the integer image displacement caused by pbc().
    1444              : !> \param snapshot ...
    1445              : !> \param cell ...
    1446              : !> \param particle_set ...
    1447              : ! **************************************************************************************************
    1448         2594 :    SUBROUTINE wfi_store_kp_pbc_shift(snapshot, cell, particle_set)
    1449              :       TYPE(qs_wf_snapshot_type), POINTER                 :: snapshot
    1450              :       TYPE(cell_type), POINTER                           :: cell
    1451              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1452              : 
    1453              :       INTEGER                                            :: iatom, natom
    1454              :       REAL(KIND=dp), DIMENSION(3)                        :: frac_pbc, frac_raw, r_pbc
    1455              : 
    1456         2594 :       CPASSERT(ASSOCIATED(snapshot))
    1457         2594 :       CPASSERT(ASSOCIATED(cell))
    1458         2594 :       CPASSERT(ASSOCIATED(particle_set))
    1459              : 
    1460         2594 :       natom = SIZE(particle_set)
    1461         2594 :       IF (ASSOCIATED(snapshot%kp_pbc_shift)) THEN
    1462          300 :          DEALLOCATE (snapshot%kp_pbc_shift)
    1463              :       END IF
    1464         7782 :       ALLOCATE (snapshot%kp_pbc_shift(3, natom))
    1465        13790 :       DO iatom = 1, natom
    1466        11196 :          r_pbc(1:3) = pbc(particle_set(iatom)%r(1:3), cell)
    1467        11196 :          CALL real_to_scaled(frac_raw, particle_set(iatom)%r(1:3), cell)
    1468        11196 :          CALL real_to_scaled(frac_pbc, r_pbc(1:3), cell)
    1469        47378 :          snapshot%kp_pbc_shift(1:3, iatom) = NINT(frac_pbc(1:3) - frac_raw(1:3))
    1470              :       END DO
    1471         2594 :    END SUBROUTINE wfi_store_kp_pbc_shift
    1472              : 
    1473              : ! **************************************************************************************************
    1474              : !> \brief Computes the current internal PBC image shift used by pbc().
    1475              : !> \param qs_env ...
    1476              : !> \param pbc_shift ...
    1477              : ! **************************************************************************************************
    1478          442 :    SUBROUTINE wfi_compute_kp_pbc_shift(qs_env, pbc_shift)
    1479              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1480              :       INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: pbc_shift
    1481              : 
    1482              :       INTEGER                                            :: iatom, natom
    1483              :       REAL(KIND=dp), DIMENSION(3)                        :: frac_pbc, frac_raw, r_pbc
    1484              :       TYPE(cell_type), POINTER                           :: cell
    1485          442 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1486              : 
    1487          442 :       NULLIFY (cell, particle_set)
    1488          442 :       CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
    1489          442 :       CPASSERT(ASSOCIATED(cell))
    1490          442 :       CPASSERT(ASSOCIATED(particle_set))
    1491              : 
    1492          442 :       natom = SIZE(particle_set)
    1493         1326 :       ALLOCATE (pbc_shift(3, natom))
    1494         3386 :       DO iatom = 1, natom
    1495         2944 :          r_pbc(1:3) = pbc(particle_set(iatom)%r(1:3), cell)
    1496         2944 :          CALL real_to_scaled(frac_raw, particle_set(iatom)%r(1:3), cell)
    1497         2944 :          CALL real_to_scaled(frac_pbc, r_pbc(1:3), cell)
    1498        12218 :          pbc_shift(1:3, iatom) = NINT(frac_pbc(1:3) - frac_raw(1:3))
    1499              :       END DO
    1500          442 :    END SUBROUTINE wfi_compute_kp_pbc_shift
    1501              : 
    1502              : ! **************************************************************************************************
    1503              : !> \brief Applies the atom-wise Bloch phase associated with a change of the internal
    1504              : !>        k-point PBC image convention to real/imaginary MO coefficient matrices.
    1505              : !> \param rmos real part of the MO coefficients
    1506              : !> \param imos imaginary part of the MO coefficients
    1507              : !> \param pbc_shift_delta target shift minus source shift for each atom
    1508              : !> \param xk fractional k-point coordinates
    1509              : !> \param matrix_template AO block structure used to map rows to atoms
    1510              : ! **************************************************************************************************
    1511         1280 :    SUBROUTINE wfi_apply_kp_pbc_phase_fm(rmos, imos, pbc_shift_delta, xk, matrix_template)
    1512              :       TYPE(cp_fm_type), POINTER                          :: rmos, imos
    1513              :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: pbc_shift_delta
    1514              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: xk
    1515              :       TYPE(dbcsr_type), POINTER                          :: matrix_template
    1516              : 
    1517              :       INTEGER                                            :: iatom, icol, irow, natom, nmo, nrow, &
    1518              :                                                             row_start
    1519         1280 :       INTEGER, DIMENSION(:), POINTER                     :: row_blk_size
    1520              :       REAL(KIND=dp)                                      :: ci, cr, i_old, r_old, theta
    1521         1280 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: iblock, rblock
    1522              : 
    1523            0 :       CPASSERT(ASSOCIATED(rmos))
    1524         1280 :       CPASSERT(ASSOCIATED(imos))
    1525         1280 :       CPASSERT(ASSOCIATED(matrix_template))
    1526              : 
    1527         1280 :       natom = SIZE(pbc_shift_delta, 2)
    1528         1280 :       CALL cp_fm_get_info(rmos, ncol_global=nmo)
    1529         1280 :       NULLIFY (row_blk_size)
    1530         1280 :       CALL dbcsr_get_info(matrix_template, row_blk_size=row_blk_size)
    1531         1280 :       CPASSERT(SIZE(row_blk_size) >= natom)
    1532              : 
    1533         1280 :       row_start = 1
    1534         7266 :       DO iatom = 1, natom
    1535         5986 :          nrow = row_blk_size(iatom)
    1536        23720 :          IF (ANY(pbc_shift_delta(1:3, iatom) /= 0)) THEN
    1537          528 :             theta = twopi*SUM(xk(1:3)*REAL(pbc_shift_delta(1:3, iatom), KIND=dp))
    1538          132 :             cr = COS(theta)
    1539          132 :             ci = SIN(theta)
    1540          792 :             ALLOCATE (rblock(nrow, nmo), iblock(nrow, nmo))
    1541          132 :             CALL cp_fm_get_submatrix(rmos, rblock, row_start, 1, nrow, nmo)
    1542          132 :             CALL cp_fm_get_submatrix(imos, iblock, row_start, 1, nrow, nmo)
    1543         4956 :             DO icol = 1, nmo
    1544        36612 :                DO irow = 1, nrow
    1545        31656 :                   r_old = rblock(irow, icol)
    1546        31656 :                   i_old = iblock(irow, icol)
    1547        31656 :                   rblock(irow, icol) = cr*r_old - ci*i_old
    1548        36480 :                   iblock(irow, icol) = ci*r_old + cr*i_old
    1549              :                END DO
    1550              :             END DO
    1551          132 :             CALL cp_fm_set_submatrix(rmos, rblock, row_start, 1, nrow, nmo)
    1552          132 :             CALL cp_fm_set_submatrix(imos, iblock, row_start, 1, nrow, nmo)
    1553          132 :             DEALLOCATE (rblock, iblock)
    1554              :          END IF
    1555         7266 :          row_start = row_start + nrow
    1556              :       END DO
    1557         2560 :    END SUBROUTINE wfi_apply_kp_pbc_phase_fm
    1558              : 
    1559              : ! **************************************************************************************************
    1560              : !> \brief Applies the atom-wise Bloch phase associated with a change of the internal
    1561              : !>        k-point PBC image convention to a complex MO coefficient matrix.
    1562              : !> \param cmos complex MO coefficients
    1563              : !> \param pbc_shift_delta target shift minus source shift for each atom
    1564              : !> \param xk fractional k-point coordinates
    1565              : !> \param matrix_template AO block structure used to map rows to atoms
    1566              : ! **************************************************************************************************
    1567         5318 :    SUBROUTINE wfi_apply_kp_pbc_phase_cfm(cmos, pbc_shift_delta, xk, matrix_template)
    1568              :       TYPE(cp_cfm_type), INTENT(INOUT)                   :: cmos
    1569              :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: pbc_shift_delta
    1570              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: xk
    1571              :       TYPE(dbcsr_type), POINTER                          :: matrix_template
    1572              : 
    1573              :       COMPLEX(KIND=dp)                                   :: phase
    1574         5318 :       COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :)     :: zblock
    1575              :       INTEGER                                            :: iatom, natom, nmo, nrow, row_start
    1576         5318 :       INTEGER, DIMENSION(:), POINTER                     :: row_blk_size
    1577              :       REAL(KIND=dp)                                      :: theta
    1578              : 
    1579            0 :       CPASSERT(ASSOCIATED(matrix_template))
    1580              : 
    1581         5318 :       natom = SIZE(pbc_shift_delta, 2)
    1582         5318 :       CALL cp_cfm_get_info(cmos, ncol_global=nmo)
    1583         5318 :       NULLIFY (row_blk_size)
    1584         5318 :       CALL dbcsr_get_info(matrix_template, row_blk_size=row_blk_size)
    1585         5318 :       CPASSERT(SIZE(row_blk_size) >= natom)
    1586              : 
    1587         5318 :       row_start = 1
    1588        37186 :       DO iatom = 1, natom
    1589        31868 :          nrow = row_blk_size(iatom)
    1590       126640 :          IF (ANY(pbc_shift_delta(1:3, iatom) /= 0)) THEN
    1591         1856 :             theta = twopi*SUM(xk(1:3)*REAL(pbc_shift_delta(1:3, iatom), KIND=dp))
    1592          464 :             phase = CMPLX(COS(theta), SIN(theta), KIND=dp)
    1593         1856 :             ALLOCATE (zblock(nrow, nmo))
    1594          464 :             CALL cp_cfm_get_submatrix(cmos, zblock, row_start, 1, nrow, nmo)
    1595       197224 :             zblock = phase*zblock
    1596          464 :             CALL cp_cfm_set_submatrix(cmos, zblock, row_start, 1, nrow, nmo)
    1597          464 :             DEALLOCATE (zblock)
    1598              :          END IF
    1599        37186 :          row_start = row_start + nrow
    1600              :       END DO
    1601        10636 :    END SUBROUTINE wfi_apply_kp_pbc_phase_cfm
    1602              : 
    1603              : ! **************************************************************************************************
    1604              : !> \brief Performs PS/ASPC wavefunction extrapolation for k-point calculations.
    1605              : !>        Applies PS/ASPC coefficients to complex MO coefficients at each k-point,
    1606              : !>        with subspace alignment via historical overlap matrices.
    1607              : !>        Delegates final orthogonalization and density building to wfi_use_prev_wf_kp.
    1608              : !> \param wf_history  wavefunction history buffer
    1609              : !> \param qs_env      QS environment
    1610              : !> \param nvec        number of history snapshots to use
    1611              : !> \param io_unit     output unit for logging
    1612              : !> \param print_level current print level
    1613              : ! **************************************************************************************************
    1614          848 :    SUBROUTINE wfi_extrapolate_ps_aspc_kp(wf_history, qs_env, nvec, io_unit, print_level)
    1615              :       TYPE(qs_wf_history_type), POINTER                  :: wf_history
    1616              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1617              :       INTEGER, INTENT(IN)                                :: nvec, io_unit, print_level
    1618              : 
    1619              :       CHARACTER(len=*), PARAMETER :: routineN = 'wfi_extrapolate_ps_aspc_kp'
    1620              : 
    1621              :       INTEGER                                            :: handle, i, ik, ikp, ispin, kplocal, &
    1622              :                                                             method_nr, nao, nmo, nspin
    1623              :       INTEGER, DIMENSION(2)                              :: kp_range
    1624              :       LOGICAL                                            :: use_real_wfn
    1625              :       REAL(KIND=dp)                                      :: alpha_coeff
    1626          424 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: xkp
    1627              :       TYPE(cp_cfm_type)                                  :: cfm_nao_nmo_work, cmos_1, cmos_i, &
    1628              :                                                             cmos_new, csc_cfm
    1629              :       TYPE(cp_fm_struct_type), POINTER                   :: nmo_nmo_struct
    1630              :       TYPE(cp_fm_type), POINTER                          :: imos, mo_coeff, rmos
    1631          424 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s_kp
    1632              :       TYPE(kpoint_env_type), POINTER                     :: kp
    1633              :       TYPE(kpoint_type), POINTER                         :: kpoints
    1634              :       TYPE(qs_wf_snapshot_type), POINTER                 :: t0_state, t1_state
    1635              : 
    1636          424 :       method_nr = wf_history%interpolation_method_nr
    1637              : 
    1638          424 :       CALL timeset(routineN, handle)
    1639          424 :       NULLIFY (kpoints, kp, mo_coeff, rmos, imos, t0_state, t1_state, nmo_nmo_struct, &
    1640          424 :                matrix_s_kp, xkp)
    1641              : 
    1642          424 :       CALL get_qs_env(qs_env, kpoints=kpoints, matrix_s_kp=matrix_s_kp)
    1643          424 :       CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn, kp_range=kp_range, xkp=xkp)
    1644          424 :       kplocal = kp_range(2) - kp_range(1) + 1
    1645              : 
    1646          424 :       IF (use_real_wfn) THEN
    1647            0 :          IF (method_nr == wfi_aspc_nr) THEN
    1648              :             CALL cp_warn(__LOCATION__, "ASPC with k-points requires complex wavefunctions; "// &
    1649            0 :                          "falling back to USE_PREV_WF.")
    1650              :          ELSE
    1651              :             CALL cp_warn(__LOCATION__, "PS with k-points requires complex wavefunctions; "// &
    1652            0 :                          "falling back to USE_PREV_WF.")
    1653              :          END IF
    1654            0 :          CALL wfi_use_prev_wf_kp(qs_env, io_unit, print_level)
    1655            0 :          CALL timestop(handle)
    1656            0 :          RETURN
    1657              :       END IF
    1658              : 
    1659          424 :       kp => kpoints%kp_env(1)%kpoint_env
    1660          424 :       nspin = SIZE(kp%mos, 2)
    1661          424 :       CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo, mo_coeff=mo_coeff)
    1662              : 
    1663          424 :       IF (method_nr == wfi_aspc_nr) THEN
    1664          420 :          IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
    1665              :             WRITE (UNIT=io_unit, FMT="(/,T2,A,/,T3,A,I0)") &
    1666           40 :                "Parameters for the always stable predictor-corrector (ASPC) method:", &
    1667           80 :                "ASPC order: ", MAX(nvec - 2, 0)
    1668              :          END IF
    1669              :       END IF
    1670              : 
    1671           16 :       IF (method_nr == wfi_aspc_nr) THEN
    1672          420 :          CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct, set_zero=.TRUE.)
    1673          420 :          CALL cp_cfm_create(cmos_1, mo_coeff%matrix_struct, set_zero=.TRUE.)
    1674          420 :          CALL cp_cfm_create(cmos_i, mo_coeff%matrix_struct, set_zero=.TRUE.)
    1675          420 :          CALL cp_cfm_create(cfm_nao_nmo_work, mo_coeff%matrix_struct, set_zero=.TRUE.)
    1676              :       ELSE
    1677            4 :          CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct)
    1678            4 :          CALL cp_cfm_create(cmos_1, mo_coeff%matrix_struct)
    1679            4 :          CALL cp_cfm_create(cmos_i, mo_coeff%matrix_struct)
    1680            4 :          CALL cp_cfm_create(cfm_nao_nmo_work, mo_coeff%matrix_struct)
    1681              :       END IF
    1682              : 
    1683              :       CALL cp_fm_struct_create(nmo_nmo_struct, template_fmstruct=mo_coeff%matrix_struct, &
    1684          424 :                                nrow_global=nmo, ncol_global=nmo)
    1685          424 :       IF (method_nr == wfi_aspc_nr) THEN
    1686          420 :          CALL cp_cfm_create(csc_cfm, nmo_nmo_struct, set_zero=.TRUE.)
    1687              :       ELSE
    1688            4 :          CALL cp_cfm_create(csc_cfm, nmo_nmo_struct)
    1689              :       END IF
    1690          424 :       CALL cp_fm_struct_release(nmo_nmo_struct)
    1691              : 
    1692              :       ! Phase 1: Initialize C_new(k) = B(1) * C_1(k)
    1693          424 :       t1_state => wfi_get_snapshot(wf_history, wf_index=1)
    1694          424 :       IF (method_nr == wfi_aspc_nr) THEN
    1695          420 :          alpha_coeff = REAL(4*nvec - 2, KIND=dp)/REAL(nvec + 1, KIND=dp)
    1696          420 :          IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
    1697           40 :             WRITE (UNIT=io_unit, FMT="(T3,A2,I0,A4,F10.6)") "B(", 1, ") = ", alpha_coeff
    1698              :          END IF
    1699              :       ELSE
    1700            4 :          alpha_coeff = nvec
    1701              :       END IF
    1702              : 
    1703         1450 :       DO ikp = 1, kplocal
    1704         1026 :          kp => kpoints%kp_env(ikp)%kpoint_env
    1705         2672 :          DO ispin = 1, nspin
    1706         1222 :             CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
    1707         1222 :             CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
    1708         1222 :             CALL cp_fm_to_fm(t1_state%wf_kp(ikp, 1, ispin), rmos)
    1709         1222 :             CALL cp_fm_to_fm(t1_state%wf_kp(ikp, 2, ispin), imos)
    1710         1222 :             CALL cp_fm_scale(alpha_coeff, rmos)
    1711         2248 :             CALL cp_fm_scale(alpha_coeff, imos)
    1712              :          END DO
    1713              :       END DO
    1714              : 
    1715              :       ! Phase 2: Accumulate historical snapshots C_new += B(i) * C_proj(k)
    1716         1736 :       DO i = 2, nvec
    1717         1312 :          t0_state => wfi_get_snapshot(wf_history, wf_index=i)
    1718         1312 :          IF (method_nr == wfi_aspc_nr) THEN
    1719              :             alpha_coeff = (-1.0_dp)**(i + 1)*REAL(i, KIND=dp)* &
    1720         1310 :                           binomial(2*nvec, nvec - i)/binomial(2*nvec - 2, nvec - 1)
    1721         1310 :             IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
    1722           71 :                WRITE (UNIT=io_unit, FMT="(T3,A2,I0,A4,F10.6)") "B(", i, ") = ", alpha_coeff
    1723              :             END IF
    1724              :          ELSE
    1725            2 :             alpha_coeff = -1.0_dp*alpha_coeff*REAL(nvec - i + 1, dp)/REAL(i, dp)
    1726              :          END IF
    1727              : 
    1728         4237 :          DO ikp = 1, kplocal
    1729         2501 :             kp => kpoints%kp_env(ikp)%kpoint_env
    1730         6410 :             DO ispin = 1, nspin
    1731         2597 :                ik = kp_range(1) + ikp - 1
    1732         2597 :                CALL cp_fm_to_cfm(t1_state%wf_kp(ikp, 1, ispin), t1_state%wf_kp(ikp, 2, ispin), cmos_1)
    1733         2597 :                CALL cp_fm_to_cfm(t0_state%wf_kp(ikp, 1, ispin), t0_state%wf_kp(ikp, 2, ispin), cmos_i)
    1734              : 
    1735              :                ! Express the reference snapshot in the image convention of snapshot i,
    1736              :                ! because the historical overlap below belongs to snapshot i.
    1737              :                CALL wfi_apply_kp_pbc_phase_cfm(cmos_1, t0_state%kp_pbc_shift - t1_state%kp_pbc_shift, &
    1738        64645 :                                                xkp(1:3, ik), matrix_s_kp(1, 1)%matrix)
    1739              : 
    1740              :                ! Subspace projection: C_proj = C_i * (C_i^dag S_i C_1)
    1741              :                CALL cp_cfm_gemm('N', 'N', nao, nmo, nao, z_one, &
    1742         2597 :                                 t0_state%overlap_cfm_kp(ikp), cmos_1, z_zero, cfm_nao_nmo_work)
    1743              :                CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, &
    1744         2597 :                                 cmos_i, cfm_nao_nmo_work, z_zero, csc_cfm)
    1745              :                CALL cp_cfm_gemm('N', 'N', nao, nmo, nmo, z_one, &
    1746         2597 :                                 cmos_i, csc_cfm, z_zero, cfm_nao_nmo_work)
    1747              : 
    1748              :                ! Convert the projected contribution from snapshot i to the reference
    1749              :                ! image convention of snapshot 1. The final conversion to the current
    1750              :                ! convention is centralized in wfi_use_prev_wf_kp.
    1751              :                CALL wfi_apply_kp_pbc_phase_cfm(cfm_nao_nmo_work, t1_state%kp_pbc_shift - t0_state%kp_pbc_shift, &
    1752        64645 :                                                xkp(1:3, ik), matrix_s_kp(1, 1)%matrix)
    1753              : 
    1754         2597 :                CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
    1755         2597 :                CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
    1756         2597 :                CALL cp_fm_to_cfm(rmos, imos, cmos_new)
    1757         2597 :                CALL cp_cfm_scale_and_add(z_one, cmos_new, CMPLX(alpha_coeff, 0.0_dp, KIND=dp), cfm_nao_nmo_work)
    1758         5098 :                CALL cp_cfm_to_fm(cmos_new, rmos, imos)
    1759              :             END DO
    1760              :          END DO
    1761              :       END DO
    1762              : 
    1763          424 :       CALL cp_cfm_release(cmos_new)
    1764          424 :       CALL cp_cfm_release(cmos_1)
    1765          424 :       CALL cp_cfm_release(cmos_i)
    1766          424 :       CALL cp_cfm_release(cfm_nao_nmo_work)
    1767          424 :       CALL cp_cfm_release(csc_cfm)
    1768              : 
    1769              :       ! Phase 3: Convert the extrapolated WFN from the reference snapshot image
    1770              :       ! convention to the current k-point PBC convention, then reorthogonalize and
    1771              :       ! rebuild the density. Keep the actual phase handling centralized in
    1772              :       ! wfi_use_prev_wf_kp so that USE_PREV_WF and ASPC/PS share the same path.
    1773              :       CALL wfi_use_prev_wf_kp(qs_env, 0, print_level, pbc_shift_ref=t1_state%kp_pbc_shift, &
    1774          424 :                               load_snapshot_wf=.FALSE.)
    1775              : 
    1776          424 :       CALL timestop(handle)
    1777              : 
    1778          424 :    END SUBROUTINE wfi_extrapolate_ps_aspc_kp
    1779              : 
    1780              : ! **************************************************************************************************
    1781              : !> \brief GEXT_PROJ/GEXT_PROJ_QTR wavefunction extrapolation for complex k-points.
    1782              : !>        This follows the existing ASPC/PS k-point projection path, but uses
    1783              : !>        the GEXT-fitted coefficients.
    1784              : !> \param wf_history wavefunction history buffer
    1785              : !> \param qs_env The QS environment
    1786              : !> \param nvec number of previous wavefunctions
    1787              : !> \param io_unit output unit
    1788              : !> \param print_level current print level
    1789              : ! **************************************************************************************************
    1790           14 :    SUBROUTINE wfi_extrapolate_gext_proj_kp(wf_history, qs_env, nvec, io_unit, print_level)
    1791              :       TYPE(qs_wf_history_type), POINTER                  :: wf_history
    1792              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1793              :       INTEGER, INTENT(IN)                                :: nvec, io_unit, print_level
    1794              : 
    1795              :       CHARACTER(len=*), PARAMETER :: routineN = 'wfi_extrapolate_gext_proj_kp'
    1796              : 
    1797              :       INTEGER                                            :: handle, i, ik, ikp, ispin, kplocal, &
    1798              :                                                             method_nr, nao, nmo, nspin
    1799              :       INTEGER, DIMENSION(2)                              :: kp_range
    1800              :       LOGICAL                                            :: use_real_wfn
    1801           14 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: coeffs, weight_kp
    1802           14 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: wkp
    1803           14 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: xkp
    1804              :       TYPE(cp_cfm_type)                                  :: cfm_nao_nmo_work, cmos_1, cmos_i, &
    1805              :                                                             cmos_new, csc_cfm
    1806           14 :       TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:)       :: csmat_cur
    1807           14 :       TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER     :: ao_ao_fm_pools_kp
    1808              :       TYPE(cp_fm_struct_type), POINTER                   :: ao_ao_struct, nmo_nmo_struct
    1809              :       TYPE(cp_fm_type), POINTER                          :: imos, mo_coeff, rmos
    1810           14 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s_kp
    1811              :       TYPE(kpoint_env_type), POINTER                     :: kp
    1812              :       TYPE(kpoint_type), POINTER                         :: kpoints
    1813              :       TYPE(mp_para_env_type), POINTER                    :: para_env_inter_kp
    1814              :       TYPE(qs_matrix_pools_type), POINTER                :: mpools_kp
    1815              :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
    1816              :       TYPE(qs_wf_snapshot_type), POINTER                 :: t0_state, t1_state
    1817              : 
    1818           14 :       method_nr = wf_history%interpolation_method_nr
    1819              : 
    1820           14 :       CALL timeset(routineN, handle)
    1821           14 :       NULLIFY (ao_ao_struct, imos, kp, kpoints, matrix_s_kp, mo_coeff, mpools_kp, &
    1822           14 :                para_env_inter_kp, rmos, scf_env, t0_state, t1_state, wkp, xkp, nmo_nmo_struct, &
    1823           14 :                ao_ao_fm_pools_kp)
    1824              : 
    1825           14 :       CALL get_qs_env(qs_env, kpoints=kpoints, matrix_s_kp=matrix_s_kp, scf_env=scf_env)
    1826              :       CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn, kp_range=kp_range, &
    1827              :                            wkp=wkp, xkp=xkp, &
    1828           14 :                            mpools=mpools_kp, para_env_inter_kp=para_env_inter_kp)
    1829           14 :       kplocal = kp_range(2) - kp_range(1) + 1
    1830              : 
    1831           14 :       IF (use_real_wfn) THEN
    1832              :          CALL cp_warn(__LOCATION__, "GExt with k-points requires complex wavefunctions; "// &
    1833            0 :                       "falling back to USE_PREV_WF.")
    1834            0 :          CALL wfi_use_prev_wf_kp(qs_env, io_unit, print_level)
    1835            0 :          CALL timestop(handle)
    1836            0 :          RETURN
    1837              :       END IF
    1838              : 
    1839           14 :       IF (nvec >= wf_history%memory_depth) THEN
    1840            0 :          IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. &
    1841              :              (qs_env%scf_control%eps_scf_hist /= 0)) THEN
    1842            0 :             qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
    1843            0 :             qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
    1844            0 :             qs_env%scf_control%outer_scf%have_scf = .FALSE.
    1845            0 :          ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
    1846            0 :             qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
    1847            0 :             qs_env%scf_control%outer_scf%have_scf = .FALSE.
    1848            0 :          ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
    1849            0 :             qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
    1850              :          END IF
    1851              :       END IF
    1852              : 
    1853           14 :       kp => kpoints%kp_env(1)%kpoint_env
    1854           14 :       nspin = SIZE(kp%mos, 2)
    1855           14 :       CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo, mo_coeff=mo_coeff)
    1856              : 
    1857           14 :       CALL mpools_get(mpools_kp, ao_ao_fm_pools=ao_ao_fm_pools_kp)
    1858           14 :       ao_ao_struct => fm_pool_get_el_struct(ao_ao_fm_pools_kp(1)%pool)
    1859              : 
    1860          108 :       ALLOCATE (csmat_cur(kplocal), weight_kp(kplocal))
    1861           52 :       DO ikp = 1, kplocal
    1862           38 :          CALL cp_cfm_create(csmat_cur(ikp), ao_ao_struct)
    1863           52 :          weight_kp(ikp) = wkp(kp_range(1) + ikp - 1)
    1864              :       END DO
    1865              : 
    1866              :       CALL kpoint_operator_get_all(kpoints, scf_env%scf_work1(1), matrix_s_kp, kpoint_slot_s, &
    1867           14 :                                    csmat_cur, matrix_row=kpoint_spin_free)
    1868              : 
    1869           42 :       ALLOCATE (coeffs(nvec))
    1870           14 :       IF (method_nr == wfi_gext_proj_nr) THEN
    1871              :          CALL diff_fitting(wf_history, matrix_s_kp(1, 1)%matrix, coeffs, nvec, &
    1872              :                            1e-4_dp, io_unit, print_level, current_overlap_kp=csmat_cur, &
    1873           10 :                            kpoint_weights=weight_kp, para_env_inter_kp=para_env_inter_kp)
    1874              :       ELSE
    1875              :          CALL tr_fitting(wf_history, matrix_s_kp(1, 1)%matrix, coeffs, nvec, &
    1876              :                          1e-4_dp, io_unit, print_level, current_overlap_kp=csmat_cur, &
    1877            4 :                          kpoint_weights=weight_kp, para_env_inter_kp=para_env_inter_kp)
    1878              :       END IF
    1879              : 
    1880              :       ! Accumulate the extrapolated WFN using the same projected-WFN path as ASPC/PS.
    1881           14 :       CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct)
    1882           14 :       CALL cp_cfm_create(cmos_1, mo_coeff%matrix_struct)
    1883           14 :       CALL cp_cfm_create(cmos_i, mo_coeff%matrix_struct)
    1884           14 :       CALL cp_cfm_create(cfm_nao_nmo_work, mo_coeff%matrix_struct)
    1885              :       CALL cp_fm_struct_create(nmo_nmo_struct, template_fmstruct=mo_coeff%matrix_struct, &
    1886           14 :                                nrow_global=nmo, ncol_global=nmo)
    1887           14 :       CALL cp_cfm_create(csc_cfm, nmo_nmo_struct)
    1888           14 :       CALL cp_fm_struct_release(nmo_nmo_struct)
    1889              : 
    1890           14 :       t1_state => wfi_get_snapshot(wf_history, wf_index=1)
    1891           52 :       DO ikp = 1, kplocal
    1892           38 :          kp => kpoints%kp_env(ikp)%kpoint_env
    1893           94 :          DO ispin = 1, nspin
    1894           42 :             CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
    1895           42 :             CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
    1896           42 :             CALL cp_fm_set_all(rmos, 0.0_dp)
    1897           80 :             CALL cp_fm_set_all(imos, 0.0_dp)
    1898              :          END DO
    1899              :       END DO
    1900              : 
    1901           34 :       DO i = 1, nvec
    1902           20 :          t0_state => wfi_get_snapshot(wf_history, wf_index=i)
    1903           90 :          DO ikp = 1, kplocal
    1904           56 :             kp => kpoints%kp_env(ikp)%kpoint_env
    1905           56 :             ik = kp_range(1) + ikp - 1
    1906          138 :             DO ispin = 1, nspin
    1907           62 :                CALL cp_fm_to_cfm(t1_state%wf_kp(ikp, 1, ispin), t1_state%wf_kp(ikp, 2, ispin), cmos_1)
    1908           62 :                CALL cp_fm_to_cfm(t0_state%wf_kp(ikp, 1, ispin), t0_state%wf_kp(ikp, 2, ispin), cmos_i)
    1909              : 
    1910              :                CALL wfi_apply_kp_pbc_phase_cfm(cmos_1, t0_state%kp_pbc_shift - t1_state%kp_pbc_shift, &
    1911         1750 :                                                xkp(1:3, ik), matrix_s_kp(1, 1)%matrix)
    1912              : 
    1913              :                CALL cp_cfm_gemm('N', 'N', nao, nmo, nao, z_one, &
    1914           62 :                                 t0_state%overlap_cfm_kp(ikp), cmos_1, z_zero, cfm_nao_nmo_work)
    1915              :                CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, &
    1916           62 :                                 cmos_i, cfm_nao_nmo_work, z_zero, csc_cfm)
    1917              :                CALL cp_cfm_gemm('N', 'N', nao, nmo, nmo, z_one, &
    1918           62 :                                 cmos_i, csc_cfm, z_zero, cfm_nao_nmo_work)
    1919              : 
    1920              :                CALL wfi_apply_kp_pbc_phase_cfm(cfm_nao_nmo_work, t1_state%kp_pbc_shift - t0_state%kp_pbc_shift, &
    1921         1750 :                                                xkp(1:3, ik), matrix_s_kp(1, 1)%matrix)
    1922              : 
    1923           62 :                CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
    1924           62 :                CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
    1925           62 :                CALL cp_fm_to_cfm(rmos, imos, cmos_new)
    1926           62 :                CALL cp_cfm_scale_and_add(z_one, cmos_new, CMPLX(coeffs(i), 0.0_dp, KIND=dp), cfm_nao_nmo_work)
    1927          118 :                CALL cp_cfm_to_fm(cmos_new, rmos, imos)
    1928              :             END DO
    1929              :          END DO
    1930              :       END DO
    1931              : 
    1932           14 :       CALL cp_cfm_release(cmos_new)
    1933           14 :       CALL cp_cfm_release(cmos_1)
    1934           14 :       CALL cp_cfm_release(cmos_i)
    1935           14 :       CALL cp_cfm_release(cfm_nao_nmo_work)
    1936           14 :       CALL cp_cfm_release(csc_cfm)
    1937              : 
    1938              :       CALL wfi_use_prev_wf_kp(qs_env, 0, print_level, pbc_shift_ref=t1_state%kp_pbc_shift, &
    1939           14 :                               load_snapshot_wf=.FALSE.)
    1940              : 
    1941           52 :       DO ikp = 1, kplocal
    1942           52 :          CALL cp_cfm_release(csmat_cur(ikp))
    1943              :       END DO
    1944           14 :       DEALLOCATE (csmat_cur, coeffs, weight_kp)
    1945              : 
    1946           14 :       CALL timestop(handle)
    1947              : 
    1948           70 :    END SUBROUTINE wfi_extrapolate_gext_proj_kp
    1949              : 
    1950              : ! **************************************************************************************************
    1951              : !> \brief Decides if scf control variables has to changed due
    1952              : !>      to using a WF extrapolation.
    1953              : !> \param qs_env The QS environment
    1954              : !> \param nvec ...
    1955              : !> \par History
    1956              : !>      11.2006 created [TdK]
    1957              : !> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch)
    1958              : ! **************************************************************************************************
    1959        11497 :    ELEMENTAL SUBROUTINE wfi_set_history_variables(qs_env, nvec)
    1960              :       TYPE(qs_environment_type), INTENT(INOUT)           :: qs_env
    1961              :       INTEGER, INTENT(IN)                                :: nvec
    1962              : 
    1963        11497 :       IF (nvec >= qs_env%wf_history%memory_depth) THEN
    1964         1897 :          IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. (qs_env%scf_control%eps_scf_hist /= 0)) THEN
    1965            0 :             qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
    1966            0 :             qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
    1967            0 :             qs_env%scf_control%outer_scf%have_scf = .FALSE.
    1968         1897 :          ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
    1969            0 :             qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
    1970            0 :             qs_env%scf_control%outer_scf%have_scf = .FALSE.
    1971         1897 :          ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
    1972            0 :             qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
    1973            0 :             qs_env%scf_control%outer_scf%eps_scf = qs_env%scf_control%eps_scf_hist
    1974              :          END IF
    1975              :       END IF
    1976              : 
    1977        11497 :    END SUBROUTINE wfi_set_history_variables
    1978              : 
    1979              : ! **************************************************************************************************
    1980              : !> \brief updates the snapshot buffer, taking a new snapshot
    1981              : !> \param wf_history the history buffer to update
    1982              : !> \param qs_env the qs_env we get the info from
    1983              : !> \param dt ...
    1984              : !> \par History
    1985              : !>      02.2003 created [fawzi]
    1986              : !> \author fawzi
    1987              : ! **************************************************************************************************
    1988        23837 :    SUBROUTINE wfi_update(wf_history, qs_env, dt)
    1989              :       TYPE(qs_wf_history_type), POINTER                  :: wf_history
    1990              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1991              :       REAL(KIND=dp), INTENT(in)                          :: dt
    1992              : 
    1993        23837 :       CPASSERT(ASSOCIATED(wf_history))
    1994        23837 :       CPASSERT(wf_history%ref_count > 0)
    1995        23837 :       CPASSERT(ASSOCIATED(qs_env))
    1996              : 
    1997        23837 :       wf_history%snapshot_count = wf_history%snapshot_count + 1
    1998        23837 :       IF (wf_history%memory_depth > 0) THEN
    1999              :          wf_history%last_state_index = MODULO(wf_history%snapshot_count, &
    2000        22616 :                                               wf_history%memory_depth) + 1
    2001              :          CALL wfs_update(snapshot=wf_history%past_states &
    2002              :                          (wf_history%last_state_index)%snapshot, wf_history=wf_history, &
    2003        22616 :                          qs_env=qs_env, dt=dt)
    2004              :       END IF
    2005        23837 :    END SUBROUTINE wfi_update
    2006              : 
    2007              : ! **************************************************************************************************
    2008              : !> \brief reorthogonalizes the mos
    2009              : !> \param qs_env the qs_env in which to orthogonalize
    2010              : !> \param v_matrix the vectors to orthogonalize
    2011              : !> \param n_col number of column of v to orthogonalize
    2012              : !> \par History
    2013              : !>      04.2003 created [fawzi]
    2014              : !> \author Fawzi Mohamed
    2015              : ! **************************************************************************************************
    2016        34518 :    SUBROUTINE reorthogonalize_vectors(qs_env, v_matrix, n_col)
    2017              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2018              :       TYPE(cp_fm_type), INTENT(IN)                       :: v_matrix
    2019              :       INTEGER, INTENT(in), OPTIONAL                      :: n_col
    2020              : 
    2021              :       CHARACTER(len=*), PARAMETER :: routineN = 'reorthogonalize_vectors'
    2022              : 
    2023              :       INTEGER                                            :: handle, my_n_col
    2024              :       LOGICAL                                            :: has_unit_metric, &
    2025              :                                                             ortho_contains_cholesky, &
    2026              :                                                             smearing_is_used
    2027              :       TYPE(cp_fm_pool_type), POINTER                     :: maxao_maxmo_fm_pool
    2028        17259 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
    2029              :       TYPE(dft_control_type), POINTER                    :: dft_control
    2030              :       TYPE(qs_matrix_pools_type), POINTER                :: mpools
    2031              :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
    2032              :       TYPE(scf_control_type), POINTER                    :: scf_control
    2033              : 
    2034        17259 :       NULLIFY (scf_env, scf_control, maxao_maxmo_fm_pool, matrix_s, mpools, dft_control)
    2035        17259 :       CALL timeset(routineN, handle)
    2036              : 
    2037        17259 :       CPASSERT(ASSOCIATED(qs_env))
    2038              : 
    2039        17259 :       CALL cp_fm_get_info(v_matrix, ncol_global=my_n_col)
    2040        17259 :       IF (PRESENT(n_col)) my_n_col = n_col
    2041              :       CALL get_qs_env(qs_env, mpools=mpools, &
    2042              :                       scf_env=scf_env, &
    2043              :                       scf_control=scf_control, &
    2044              :                       matrix_s=matrix_s, &
    2045        17259 :                       dft_control=dft_control)
    2046        17259 :       CALL mpools_get(mpools, maxao_maxmo_fm_pool=maxao_maxmo_fm_pool)
    2047        17259 :       IF (ASSOCIATED(scf_env)) THEN
    2048              :          ortho_contains_cholesky = (scf_env%method /= ot_method_nr) .AND. &
    2049              :                                    (scf_env%cholesky_method > 0) .AND. &
    2050        17259 :                                    ASSOCIATED(scf_env%ortho)
    2051              :       ELSE
    2052              :          ortho_contains_cholesky = .FALSE.
    2053              :       END IF
    2054              : 
    2055        17259 :       CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
    2056        17259 :       smearing_is_used = .FALSE.
    2057        17259 :       IF (dft_control%smear) THEN
    2058         1906 :          smearing_is_used = .TRUE.
    2059              :       END IF
    2060              : 
    2061        17259 :       IF (has_unit_metric) THEN
    2062         3410 :          CALL make_basis_simple(v_matrix, my_n_col)
    2063        13849 :       ELSE IF (smearing_is_used) THEN
    2064              :          CALL make_basis_lowdin(vmatrix=v_matrix, ncol=my_n_col, &
    2065         1906 :                                 matrix_s=matrix_s(1)%matrix)
    2066        11943 :       ELSE IF (ortho_contains_cholesky) THEN
    2067              :          CALL make_basis_cholesky(vmatrix=v_matrix, ncol=my_n_col, &
    2068         8234 :                                   ortho=scf_env%ortho)
    2069              :       ELSE
    2070         3709 :          CALL make_basis_sm(v_matrix, my_n_col, matrix_s(1)%matrix)
    2071              :       END IF
    2072        17259 :       CALL timestop(handle)
    2073        17259 :    END SUBROUTINE reorthogonalize_vectors
    2074              : 
    2075              : ! **************************************************************************************************
    2076              : !> \brief purges wf_history retaining only the latest snapshot
    2077              : !> \param qs_env the qs env with the latest result, and that will contain
    2078              : !>        the purged wf_history
    2079              : !> \par History
    2080              : !>      05.2016 created [Nico Holmberg]
    2081              : !> \author Nico Holmberg
    2082              : ! **************************************************************************************************
    2083            0 :    SUBROUTINE wfi_purge_history(qs_env)
    2084              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2085              : 
    2086              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'wfi_purge_history'
    2087              : 
    2088              :       INTEGER                                            :: handle, io_unit, print_level
    2089              :       TYPE(cp_logger_type), POINTER                      :: logger
    2090              :       TYPE(dft_control_type), POINTER                    :: dft_control
    2091              :       TYPE(qs_wf_history_type), POINTER                  :: wf_history
    2092              : 
    2093            0 :       NULLIFY (dft_control, wf_history)
    2094              : 
    2095            0 :       CALL timeset(routineN, handle)
    2096            0 :       logger => cp_get_default_logger()
    2097            0 :       print_level = logger%iter_info%print_level
    2098              :       io_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%SCF%PRINT%PROGRAM_RUN_INFO", &
    2099            0 :                                      extension=".scfLog")
    2100              : 
    2101            0 :       CPASSERT(ASSOCIATED(qs_env))
    2102            0 :       CPASSERT(ASSOCIATED(qs_env%wf_history))
    2103            0 :       CPASSERT(qs_env%wf_history%ref_count > 0)
    2104            0 :       CALL get_qs_env(qs_env, dft_control=dft_control)
    2105              : 
    2106            0 :       SELECT CASE (qs_env%wf_history%interpolation_method_nr)
    2107              :       CASE (wfi_use_guess_method_nr, wfi_use_prev_wf_method_nr, &
    2108              :             wfi_use_prev_p_method_nr, wfi_frozen_method_nr)
    2109              :          ! do nothing
    2110              :       CASE (wfi_linear_wf_method_nr, wfi_linear_p_method_nr, &
    2111              :             wfi_linear_ps_method_nr, wfi_ps_method_nr, &
    2112              :             wfi_aspc_nr, wfi_gext_proj_nr, wfi_gext_proj_qtr_nr)
    2113            0 :          IF (qs_env%wf_history%snapshot_count >= 2) THEN
    2114            0 :             IF (debug_this_module .AND. io_unit > 0) THEN
    2115            0 :                WRITE (io_unit, FMT="(T2,A)") "QS| Purging WFN history"
    2116              :             END IF
    2117              :             CALL wfi_create(wf_history, interpolation_method_nr= &
    2118              :                             dft_control%qs_control%wf_interpolation_method_nr, &
    2119              :                             extrapolation_order=dft_control%qs_control%wf_extrapolation_order, &
    2120            0 :                             has_unit_metric=qs_env%has_unit_metric)
    2121              :             CALL set_qs_env(qs_env=qs_env, &
    2122            0 :                             wf_history=wf_history)
    2123            0 :             CALL wfi_release(wf_history)
    2124            0 :             CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
    2125              :          END IF
    2126              :       CASE DEFAULT
    2127            0 :          CPABORT("Unknown extrapolation method.")
    2128              :       END SELECT
    2129            0 :       CALL timestop(handle)
    2130              : 
    2131            0 :    END SUBROUTINE wfi_purge_history
    2132              : 
    2133              : ! **************************************************************************************************
    2134              : !> \brief Gives the coefficients that best approximate the new overlap
    2135              : !>        as a linear combination of the previous overlaps in the
    2136              : !>        wf_history buffer. This is done by solving
    2137              : !>        argmin_a || S_{n+1} - S_{n} - \sum_i^{nvec-1} a_i (S_{n-q+i} - S_{n}) ||^2
    2138              : !> \param wf_history wavefunction history buffer, containing the previous overlaps
    2139              : !> \param current_overlap current overlap in dbcsr format
    2140              : !> \param coeffs resulting nvec coefficients
    2141              : !> \param nvec number of previous overlaps
    2142              : !> \param eps Tikhonov regularization
    2143              : !> \param io_unit output unit
    2144              : !> \param print_level print level
    2145              : !> \param current_overlap_kp ...
    2146              : !> \param kpoint_weights ...
    2147              : !> \param para_env_inter_kp ...
    2148              : !> \par History
    2149              : !>      04.2026 created [Michele Nottoli]
    2150              : !> \author Michele Nottoli
    2151              : ! **************************************************************************************************
    2152          146 :    SUBROUTINE diff_fitting(wf_history, current_overlap, coeffs, nvec, eps, io_unit, print_level, &
    2153          146 :                            current_overlap_kp, kpoint_weights, para_env_inter_kp)
    2154              :       TYPE(qs_wf_history_type), POINTER                  :: wf_history
    2155              :       TYPE(dbcsr_type), INTENT(IN)                       :: current_overlap
    2156              :       INTEGER, INTENT(IN)                                :: nvec
    2157              :       REAL(KIND=dp), INTENT(OUT)                         :: coeffs(nvec)
    2158              :       REAL(KIND=dp), INTENT(IN)                          :: eps
    2159              :       INTEGER, INTENT(IN)                                :: io_unit, print_level
    2160              :       TYPE(cp_cfm_type), DIMENSION(:), INTENT(IN), &
    2161              :          OPTIONAL                                        :: current_overlap_kp
    2162              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL  :: kpoint_weights
    2163              :       TYPE(mp_para_env_type), OPTIONAL, POINTER          :: para_env_inter_kp
    2164              : 
    2165              :       COMPLEX(KIND=dp)                                   :: ztrace
    2166              :       INTEGER                                            :: i, icol_local, ikp, info, irow_local, j
    2167              :       REAL(KIND=dp)                                      :: error, norm_ref, weight
    2168          146 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: b
    2169          146 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: A
    2170              :       TYPE(cp_cfm_type)                                  :: target_diff_cfm, tmp_conj_cfm, &
    2171              :                                                             tmp_i_cfm, tmp_j_cfm
    2172              :       TYPE(dbcsr_type)                                   :: target_diff, tmp_i, tmp_j, tmp_k
    2173              :       TYPE(qs_wf_snapshot_type), POINTER                 :: ref_state, state
    2174              : 
    2175          146 :       IF (nvec <= 0) THEN
    2176            0 :          CPABORT("Not enough vectors to do the fitting")
    2177          146 :       ELSE IF (nvec == 1) THEN
    2178           26 :          coeffs(1) = 1.0_dp
    2179           86 :          RETURN
    2180              :       END IF
    2181              : 
    2182          120 :       IF (PRESENT(current_overlap_kp)) THEN
    2183           24 :          ALLOCATE (A(nvec - 1, nvec - 1), b(nvec - 1))
    2184            4 :          A = 0.0_dp
    2185            4 :          b = 0.0_dp
    2186              : 
    2187            4 :          ref_state => wfi_get_snapshot(wf_history, wf_index=1)
    2188            4 :          CALL cp_cfm_create(target_diff_cfm, current_overlap_kp(1)%matrix_struct)
    2189            4 :          CALL cp_cfm_create(tmp_i_cfm, current_overlap_kp(1)%matrix_struct)
    2190            4 :          CALL cp_cfm_create(tmp_j_cfm, current_overlap_kp(1)%matrix_struct)
    2191            4 :          CALL cp_cfm_create(tmp_conj_cfm, current_overlap_kp(1)%matrix_struct)
    2192              : 
    2193           14 :          DO ikp = 1, SIZE(current_overlap_kp)
    2194           10 :             weight = 1.0_dp
    2195           10 :             IF (PRESENT(kpoint_weights)) weight = kpoint_weights(ikp)
    2196              : 
    2197           10 :             CALL cp_cfm_to_cfm(current_overlap_kp(ikp), target_diff_cfm)
    2198              :             CALL cp_cfm_scale_and_add(z_one, target_diff_cfm, CMPLX(-1.0_dp, 0.0_dp, KIND=dp), &
    2199           10 :                                       ref_state%overlap_cfm_kp(ikp))
    2200           24 :             DO i = 2, nvec
    2201           10 :                state => wfi_get_snapshot(wf_history, wf_index=i)
    2202           10 :                CALL cp_cfm_to_cfm(state%overlap_cfm_kp(ikp), tmp_i_cfm)
    2203              :                CALL cp_cfm_scale_and_add(z_one, tmp_i_cfm, CMPLX(-1.0_dp, 0.0_dp, KIND=dp), &
    2204           10 :                                          ref_state%overlap_cfm_kp(ikp))
    2205           10 :                CALL cp_cfm_to_cfm(tmp_i_cfm, tmp_conj_cfm)
    2206          272 :                DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
    2207         4386 :                   DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
    2208              :                      tmp_conj_cfm%local_data(irow_local, icol_local) = &
    2209         4376 :                         CONJG(tmp_conj_cfm%local_data(irow_local, icol_local))
    2210              :                   END DO
    2211              :                END DO
    2212           10 :                CALL cp_cfm_trace(tmp_conj_cfm, target_diff_cfm, ztrace)
    2213           10 :                b(i - 1) = b(i - 1) + weight*REAL(ztrace, KIND=dp)
    2214              : 
    2215           30 :                DO j = 2, i
    2216           10 :                   state => wfi_get_snapshot(wf_history, wf_index=j)
    2217           10 :                   CALL cp_cfm_to_cfm(state%overlap_cfm_kp(ikp), tmp_j_cfm)
    2218              :                   CALL cp_cfm_scale_and_add(z_one, tmp_j_cfm, CMPLX(-1.0_dp, 0.0_dp, KIND=dp), &
    2219           10 :                                             ref_state%overlap_cfm_kp(ikp))
    2220           10 :                   CALL cp_cfm_trace(tmp_conj_cfm, tmp_j_cfm, ztrace)
    2221           20 :                   A(j - 1, i - 1) = A(j - 1, i - 1) + weight*REAL(ztrace, KIND=dp)
    2222              :                END DO
    2223              :             END DO
    2224              :          END DO
    2225              : 
    2226            8 :          DO i = 2, nvec
    2227           12 :             DO j = 2, i
    2228            8 :                A(i - 1, j - 1) = A(j - 1, i - 1)
    2229              :             END DO
    2230              :          END DO
    2231              : 
    2232            4 :          IF (PRESENT(para_env_inter_kp)) THEN
    2233            4 :             IF (ASSOCIATED(para_env_inter_kp)) THEN
    2234            4 :                CALL para_env_inter_kp%sum(A)
    2235            4 :                CALL para_env_inter_kp%sum(b)
    2236              :             END IF
    2237              :          END IF
    2238              : 
    2239            8 :          DO i = 1, nvec - 1
    2240            8 :             A(i, i) = A(i, i) + eps**2
    2241              :          END DO
    2242              : 
    2243            4 :          CALL dposv('u', nvec - 1, 1, A, nvec - 1, b, nvec - 1, info)
    2244            4 :          IF (info /= 0) THEN
    2245            0 :             CPABORT("DPOSV failed.")
    2246              :          END IF
    2247              : 
    2248            8 :          coeffs(1) = 1.0_dp - SUM(b)
    2249            8 :          coeffs(2:nvec) = b(:)
    2250              : 
    2251            4 :          IF (print_level > low_print_level) THEN
    2252            2 :             error = 0.0_dp
    2253            2 :             norm_ref = 0.0_dp
    2254           10 :             DO ikp = 1, SIZE(current_overlap_kp)
    2255            8 :                weight = 1.0_dp
    2256            8 :                IF (PRESENT(kpoint_weights)) weight = kpoint_weights(ikp)
    2257            8 :                CALL cp_cfm_to_cfm(current_overlap_kp(ikp), tmp_i_cfm)
    2258           24 :                DO i = 1, nvec
    2259           16 :                   state => wfi_get_snapshot(wf_history, wf_index=i)
    2260              :                   CALL cp_cfm_scale_and_add(z_one, tmp_i_cfm, CMPLX(-coeffs(i), 0.0_dp, KIND=dp), &
    2261           24 :                                             state%overlap_cfm_kp(ikp))
    2262              :                END DO
    2263            8 :                CALL cp_cfm_to_cfm(tmp_i_cfm, tmp_conj_cfm)
    2264          264 :                DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
    2265         4360 :                   DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
    2266              :                      tmp_conj_cfm%local_data(irow_local, icol_local) = &
    2267         4352 :                         CONJG(tmp_conj_cfm%local_data(irow_local, icol_local))
    2268              :                   END DO
    2269              :                END DO
    2270            8 :                CALL cp_cfm_trace(tmp_conj_cfm, tmp_i_cfm, ztrace)
    2271            8 :                error = error + weight*REAL(ztrace, KIND=dp)
    2272            8 :                CALL cp_cfm_to_cfm(ref_state%overlap_cfm_kp(ikp), tmp_conj_cfm)
    2273          264 :                DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
    2274         4360 :                   DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
    2275              :                      tmp_conj_cfm%local_data(irow_local, icol_local) = &
    2276         4352 :                         CONJG(tmp_conj_cfm%local_data(irow_local, icol_local))
    2277              :                   END DO
    2278              :                END DO
    2279            8 :                CALL cp_cfm_trace(tmp_conj_cfm, ref_state%overlap_cfm_kp(ikp), ztrace)
    2280           18 :                norm_ref = norm_ref + weight*REAL(ztrace, KIND=dp)
    2281              :             END DO
    2282            2 :             IF (PRESENT(para_env_inter_kp)) THEN
    2283            2 :                IF (ASSOCIATED(para_env_inter_kp)) THEN
    2284            2 :                   CALL para_env_inter_kp%sum(error)
    2285            2 :                   CALL para_env_inter_kp%sum(norm_ref)
    2286              :                END IF
    2287              :             END IF
    2288            2 :             IF (io_unit > 0) THEN
    2289            1 :                WRITE (UNIT=io_unit, FMT="(/,T2,A,F20.10)") "GEXT overlap fitting error:", &
    2290            2 :                   SQRT(error/MAX(norm_ref, TINY(1.0_dp)))
    2291              :             END IF
    2292              :          END IF
    2293              : 
    2294            4 :          CALL cp_cfm_release(target_diff_cfm)
    2295            4 :          CALL cp_cfm_release(tmp_i_cfm)
    2296            4 :          CALL cp_cfm_release(tmp_j_cfm)
    2297            4 :          CALL cp_cfm_release(tmp_conj_cfm)
    2298            4 :          DEALLOCATE (A, b)
    2299            4 :          RETURN
    2300              :       END IF
    2301              : 
    2302          696 :       ALLOCATE (A(nvec - 1, nvec - 1), b(nvec - 1))
    2303              : 
    2304              :       ! get the reference for the difference fitting
    2305          116 :       ref_state => wfi_get_snapshot(wf_history, wf_index=1)
    2306              : 
    2307              :       ! assemble the target difference
    2308          116 :       CALL dbcsr_copy(target_diff, current_overlap)
    2309          116 :       CALL dbcsr_add(target_diff, ref_state%overlap, 1.0_dp, -1.0_dp)
    2310              : 
    2311              :       ! allocate tmp_k
    2312          116 :       CALL dbcsr_copy(tmp_k, current_overlap)
    2313              : 
    2314              :       ! assemble the matrix A and the RHS b
    2315          348 :       DO i = 2, nvec
    2316          232 :          state => wfi_get_snapshot(wf_history, wf_index=i)
    2317          232 :          CALL dbcsr_copy(tmp_i, state%overlap)
    2318          232 :          CALL dbcsr_add(tmp_i, ref_state%overlap, 1.0_dp, -1.0_dp)
    2319          232 :          CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_i, target_diff, 0.0_dp, tmp_k)
    2320          232 :          CALL dbcsr_trace(tmp_k, b(i - 1))
    2321              : 
    2322          724 :          DO j = 2, i
    2323          376 :             state => wfi_get_snapshot(wf_history, wf_index=j)
    2324          376 :             CALL dbcsr_copy(tmp_j, state%overlap)
    2325          376 :             CALL dbcsr_add(tmp_j, ref_state%overlap, 1.0_dp, -1.0_dp)
    2326          376 :             CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_i, tmp_j, 0.0_dp, tmp_k)
    2327          376 :             CALL dbcsr_trace(tmp_k, A(j - 1, i - 1))
    2328          608 :             A(i - 1, j - 1) = A(j - 1, i - 1)
    2329              :          END DO
    2330              :       END DO
    2331              : 
    2332              :       ! add the Tikhonov regularization
    2333          348 :       DO i = 1, nvec - 1
    2334          348 :          A(i, i) = A(i, i) + eps**2
    2335              :       END DO
    2336              : 
    2337              :       ! solve the linear system
    2338          116 :       CALL dposv('u', nvec - 1, 1, A, nvec - 1, b, nvec - 1, info)
    2339          116 :       IF (info /= 0) THEN
    2340            0 :          CPABORT("DPOSV failed.")
    2341              :       END IF
    2342              : 
    2343              :       ! set the coefficient for the reference snapshot
    2344          348 :       coeffs(1) = 1.0_dp - SUM(b)
    2345          348 :       coeffs(2:nvec) = b(:)
    2346              : 
    2347              :       ! as a consistency check, print how well the current overlap
    2348              :       ! is approximated by the linear combination of previous overlaps
    2349          116 :       IF (print_level > low_print_level) THEN
    2350           20 :          CALL dbcsr_copy(tmp_i, current_overlap)
    2351           96 :          DO i = 1, nvec
    2352           76 :             state => wfi_get_snapshot(wf_history, wf_index=i)
    2353           96 :             CALL dbcsr_add(tmp_i, state%overlap, 1.0_dp, -coeffs(i))
    2354              :          END DO
    2355           20 :          error = dbcsr_frobenius_norm(tmp_i)/dbcsr_frobenius_norm(state%overlap)
    2356           20 :          IF (io_unit > 0) THEN
    2357           10 :             WRITE (UNIT=io_unit, FMT="(/,T2,A,F20.10)") "GEXT overlap fitting error:", error
    2358              :          END IF
    2359              :       END IF
    2360              : 
    2361              :       ! free the memory
    2362          116 :       CALL dbcsr_release(tmp_i)
    2363          116 :       CALL dbcsr_release(tmp_j)
    2364          116 :       CALL dbcsr_release(tmp_k)
    2365          116 :       CALL dbcsr_release(target_diff)
    2366          116 :       DEALLOCATE (A, b)
    2367              : 
    2368          206 :    END SUBROUTINE diff_fitting
    2369              : 
    2370              : ! **************************************************************************************************
    2371              : !> \brief Gives the coefficients that best approximate the new overlap
    2372              : !>        as a time reversible linear combination of the previous overlaps in the
    2373              : !>        wf_history buffer. This is done by solving
    2374              : !>        argmin_a || S_{n+1} + S_{n+1-nvec}
    2375              : !>                    - \sum_{i=1}^q a_i (S_{n+1-nvec+i} + S_{n+1-i}) ||^2
    2376              : !>        with q = nvec/2 if nvec is even, or q = (nvec-1)/2 if odd.
    2377              : !> \param wf_history wavefunction history buffer, containing the previous overlaps
    2378              : !> \param current_overlap current overlap in dbcsr format
    2379              : !> \param coeffs resulting nvec coefficients
    2380              : !> \param nvec number of previous overlaps
    2381              : !> \param eps Tikhonov regularization
    2382              : !> \param io_unit output unit
    2383              : !> \param print_level print level
    2384              : !> \param current_overlap_kp ...
    2385              : !> \param kpoint_weights ...
    2386              : !> \param para_env_inter_kp ...
    2387              : !> \par History
    2388              : !>      04.2026 created [Michele Nottoli]
    2389              : ! **************************************************************************************************
    2390           28 :    SUBROUTINE tr_fitting(wf_history, current_overlap, coeffs, nvec, eps, io_unit, print_level, &
    2391           28 :                          current_overlap_kp, kpoint_weights, para_env_inter_kp)
    2392              :       TYPE(qs_wf_history_type), POINTER                  :: wf_history
    2393              :       TYPE(dbcsr_type), INTENT(IN)                       :: current_overlap
    2394              :       INTEGER, INTENT(IN)                                :: nvec
    2395              :       REAL(KIND=dp), INTENT(OUT)                         :: coeffs(nvec)
    2396              :       REAL(KIND=dp), INTENT(IN)                          :: eps
    2397              :       INTEGER, INTENT(IN)                                :: io_unit, print_level
    2398              :       TYPE(cp_cfm_type), DIMENSION(:), INTENT(IN), &
    2399              :          OPTIONAL                                        :: current_overlap_kp
    2400              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL  :: kpoint_weights
    2401              :       TYPE(mp_para_env_type), OPTIONAL, POINTER          :: para_env_inter_kp
    2402              : 
    2403              :       COMPLEX(KIND=dp)                                   :: ztrace
    2404              :       INTEGER                                            :: i, icol_local, ikp, info, irow_local, j, &
    2405              :                                                             ntr
    2406              :       REAL(KIND=dp)                                      :: error, norm_ref, weight
    2407           28 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: b
    2408           28 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: A
    2409              :       TYPE(cp_cfm_type)                                  :: target_overlap_cfm, tmp_conj_cfm, &
    2410              :                                                             tmp_i_cfm, tmp_j_cfm
    2411              :       TYPE(dbcsr_type)                                   :: target_overlap, tmp_i, tmp_j, tmp_k
    2412              :       TYPE(qs_wf_snapshot_type), POINTER                 :: ref_state, state
    2413              : 
    2414           28 :       IF (nvec <= 0) THEN
    2415            0 :          CPABORT("Not enough vectors to do the fitting")
    2416           28 :       ELSE IF (nvec == 1) THEN
    2417            6 :          coeffs(1) = 1.0_dp
    2418           22 :          RETURN
    2419              :       END IF
    2420              : 
    2421           22 :       IF (MOD(nvec, 2) == 0) THEN
    2422           10 :          ntr = nvec/2
    2423              :       ELSE
    2424           12 :          ntr = (nvec - 1)/2
    2425              :       END IF
    2426              : 
    2427           22 :       IF (PRESENT(current_overlap_kp)) THEN
    2428           12 :          ALLOCATE (A(ntr, ntr), b(ntr))
    2429            2 :          A = 0.0_dp
    2430            2 :          b = 0.0_dp
    2431              : 
    2432            2 :          ref_state => wfi_get_snapshot(wf_history, wf_index=nvec)
    2433            2 :          CALL cp_cfm_create(target_overlap_cfm, current_overlap_kp(1)%matrix_struct)
    2434            2 :          CALL cp_cfm_create(tmp_i_cfm, current_overlap_kp(1)%matrix_struct)
    2435            2 :          CALL cp_cfm_create(tmp_j_cfm, current_overlap_kp(1)%matrix_struct)
    2436            2 :          CALL cp_cfm_create(tmp_conj_cfm, current_overlap_kp(1)%matrix_struct)
    2437              : 
    2438           10 :          DO ikp = 1, SIZE(current_overlap_kp)
    2439            8 :             weight = 1.0_dp
    2440            8 :             IF (PRESENT(kpoint_weights)) weight = kpoint_weights(ikp)
    2441              : 
    2442            8 :             CALL cp_cfm_to_cfm(current_overlap_kp(ikp), target_overlap_cfm)
    2443            8 :             CALL cp_cfm_scale_and_add(z_one, target_overlap_cfm, z_one, ref_state%overlap_cfm_kp(ikp))
    2444           18 :             DO i = 1, ntr
    2445            8 :                state => wfi_get_snapshot(wf_history, wf_index=i)
    2446            8 :                CALL cp_cfm_to_cfm(state%overlap_cfm_kp(ikp), tmp_i_cfm)
    2447            8 :                state => wfi_get_snapshot(wf_history, wf_index=nvec - i)
    2448            8 :                CALL cp_cfm_scale_and_add(z_one, tmp_i_cfm, z_one, state%overlap_cfm_kp(ikp))
    2449              : 
    2450            8 :                CALL cp_cfm_to_cfm(tmp_i_cfm, tmp_conj_cfm)
    2451          264 :                DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
    2452         4360 :                   DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
    2453              :                      tmp_conj_cfm%local_data(irow_local, icol_local) = &
    2454         4352 :                         CONJG(tmp_conj_cfm%local_data(irow_local, icol_local))
    2455              :                   END DO
    2456              :                END DO
    2457            8 :                CALL cp_cfm_trace(tmp_conj_cfm, target_overlap_cfm, ztrace)
    2458            8 :                b(i) = b(i) + weight*REAL(ztrace, KIND=dp)
    2459           24 :                DO j = 1, i
    2460            8 :                   state => wfi_get_snapshot(wf_history, wf_index=j)
    2461            8 :                   CALL cp_cfm_to_cfm(state%overlap_cfm_kp(ikp), tmp_j_cfm)
    2462            8 :                   state => wfi_get_snapshot(wf_history, wf_index=nvec - j)
    2463            8 :                   CALL cp_cfm_scale_and_add(z_one, tmp_j_cfm, z_one, state%overlap_cfm_kp(ikp))
    2464            8 :                   CALL cp_cfm_trace(tmp_conj_cfm, tmp_j_cfm, ztrace)
    2465           16 :                   A(j, i) = A(j, i) + weight*REAL(ztrace, KIND=dp)
    2466              :                END DO
    2467              :             END DO
    2468              :          END DO
    2469              : 
    2470            4 :          DO i = 1, ntr
    2471            6 :             DO j = 1, i
    2472            4 :                A(i, j) = A(j, i)
    2473              :             END DO
    2474              :          END DO
    2475              : 
    2476            2 :          IF (PRESENT(para_env_inter_kp)) THEN
    2477            2 :             IF (ASSOCIATED(para_env_inter_kp)) THEN
    2478            2 :                CALL para_env_inter_kp%sum(A)
    2479            2 :                CALL para_env_inter_kp%sum(b)
    2480              :             END IF
    2481              :          END IF
    2482              : 
    2483            4 :          DO i = 1, ntr
    2484            4 :             A(i, i) = A(i, i) + eps**2
    2485              :          END DO
    2486              : 
    2487            2 :          CALL dposv('u', ntr, 1, A, ntr, b, ntr, info)
    2488            2 :          IF (info /= 0) THEN
    2489            0 :             CPABORT("DPOSV failed.")
    2490              :          END IF
    2491              : 
    2492            6 :          coeffs = 0.0_dp
    2493            2 :          coeffs(nvec) = -1.0_dp
    2494            4 :          DO i = 1, ntr
    2495            2 :             coeffs(i) = coeffs(i) + b(i)
    2496            4 :             coeffs(nvec - i) = coeffs(nvec - i) + b(i)
    2497              :          END DO
    2498              : 
    2499            2 :          IF (print_level > low_print_level) THEN
    2500            2 :             error = 0.0_dp
    2501            2 :             norm_ref = 0.0_dp
    2502           10 :             DO ikp = 1, SIZE(current_overlap_kp)
    2503            8 :                weight = 1.0_dp
    2504            8 :                IF (PRESENT(kpoint_weights)) weight = kpoint_weights(ikp)
    2505            8 :                CALL cp_cfm_to_cfm(current_overlap_kp(ikp), tmp_i_cfm)
    2506           24 :                DO i = 1, nvec
    2507           16 :                   state => wfi_get_snapshot(wf_history, wf_index=i)
    2508              :                   CALL cp_cfm_scale_and_add(z_one, tmp_i_cfm, CMPLX(-coeffs(i), 0.0_dp, KIND=dp), &
    2509           24 :                                             state%overlap_cfm_kp(ikp))
    2510              :                END DO
    2511            8 :                CALL cp_cfm_to_cfm(tmp_i_cfm, tmp_conj_cfm)
    2512          264 :                DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
    2513         4360 :                   DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
    2514              :                      tmp_conj_cfm%local_data(irow_local, icol_local) = &
    2515         4352 :                         CONJG(tmp_conj_cfm%local_data(irow_local, icol_local))
    2516              :                   END DO
    2517              :                END DO
    2518            8 :                CALL cp_cfm_trace(tmp_conj_cfm, tmp_i_cfm, ztrace)
    2519            8 :                error = error + weight*REAL(ztrace, KIND=dp)
    2520            8 :                CALL cp_cfm_to_cfm(ref_state%overlap_cfm_kp(ikp), tmp_conj_cfm)
    2521          264 :                DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
    2522         4360 :                   DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
    2523              :                      tmp_conj_cfm%local_data(irow_local, icol_local) = &
    2524         4352 :                         CONJG(tmp_conj_cfm%local_data(irow_local, icol_local))
    2525              :                   END DO
    2526              :                END DO
    2527            8 :                CALL cp_cfm_trace(tmp_conj_cfm, ref_state%overlap_cfm_kp(ikp), ztrace)
    2528           18 :                norm_ref = norm_ref + weight*REAL(ztrace, KIND=dp)
    2529              :             END DO
    2530            2 :             IF (PRESENT(para_env_inter_kp)) THEN
    2531            2 :                IF (ASSOCIATED(para_env_inter_kp)) THEN
    2532            2 :                   CALL para_env_inter_kp%sum(error)
    2533            2 :                   CALL para_env_inter_kp%sum(norm_ref)
    2534              :                END IF
    2535              :             END IF
    2536            2 :             IF (io_unit > 0) THEN
    2537            1 :                WRITE (UNIT=io_unit, FMT="(/,T2,A,F20.10)") "GEXT overlap fitting error:", &
    2538            2 :                   SQRT(error/MAX(norm_ref, TINY(1.0_dp)))
    2539              :             END IF
    2540              :          END IF
    2541              : 
    2542            2 :          CALL cp_cfm_release(target_overlap_cfm)
    2543            2 :          CALL cp_cfm_release(tmp_i_cfm)
    2544            2 :          CALL cp_cfm_release(tmp_j_cfm)
    2545            2 :          CALL cp_cfm_release(tmp_conj_cfm)
    2546            2 :          DEALLOCATE (A, b)
    2547            2 :          RETURN
    2548              :       END IF
    2549              : 
    2550          120 :       ALLOCATE (A(ntr, ntr), b(ntr))
    2551              : 
    2552              :       ! get the reference for the difference fitting
    2553           20 :       ref_state => wfi_get_snapshot(wf_history, wf_index=nvec)
    2554              : 
    2555              :       ! assemble the target sum
    2556           20 :       CALL dbcsr_copy(target_overlap, current_overlap)
    2557           20 :       CALL dbcsr_add(target_overlap, ref_state%overlap, 1.0_dp, 1.0_dp)
    2558              : 
    2559              :       ! allocate tmp_k
    2560           20 :       CALL dbcsr_copy(tmp_k, current_overlap)
    2561              : 
    2562              :       ! assemble the matrix A and the RHS b
    2563           52 :       DO i = 1, ntr
    2564           32 :          state => wfi_get_snapshot(wf_history, wf_index=i)
    2565           32 :          CALL dbcsr_copy(tmp_i, state%overlap)
    2566           32 :          state => wfi_get_snapshot(wf_history, wf_index=nvec - i)
    2567           32 :          CALL dbcsr_add(tmp_i, state%overlap, 1.0_dp, 1.0_dp)
    2568              : 
    2569              :          CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_i, target_overlap, &
    2570           32 :                              0.0_dp, tmp_k)
    2571           32 :          CALL dbcsr_trace(tmp_k, b(i))
    2572           96 :          DO j = 1, i
    2573           44 :             state => wfi_get_snapshot(wf_history, wf_index=j)
    2574           44 :             CALL dbcsr_copy(tmp_j, state%overlap)
    2575           44 :             state => wfi_get_snapshot(wf_history, wf_index=nvec - j)
    2576           44 :             CALL dbcsr_add(tmp_j, state%overlap, 1.0_dp, 1.0_dp)
    2577           44 :             CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_i, tmp_j, 0.0_dp, tmp_k)
    2578           44 :             CALL dbcsr_trace(tmp_k, A(j, i))
    2579           76 :             A(i, j) = A(j, i)
    2580              :          END DO
    2581              :       END DO
    2582              : 
    2583              :       ! add the Tikhonov regularization
    2584           52 :       DO i = 1, ntr
    2585           52 :          A(i, i) = A(i, i) + eps**2
    2586              :       END DO
    2587              : 
    2588              :       ! solve the linear system
    2589           20 :       CALL dposv('u', ntr, 1, A, ntr, b, ntr, info)
    2590           20 :       IF (info /= 0) THEN
    2591            0 :          CPABORT("DPOSV failed.")
    2592              :       END IF
    2593              : 
    2594              :       ! reorder the coefficients
    2595           96 :       coeffs = 0.0_dp
    2596           20 :       coeffs(nvec) = -1.0_dp
    2597           52 :       DO i = 1, ntr
    2598           32 :          coeffs(i) = coeffs(i) + b(i)
    2599           52 :          coeffs(nvec - i) = coeffs(nvec - i) + b(i)
    2600              :       END DO
    2601              : 
    2602              :       ! as a consistency check, print how well the current overlap
    2603              :       ! is approximated by the linear combination of previous overlaps
    2604           20 :       IF (print_level > low_print_level) THEN
    2605           20 :          CALL dbcsr_copy(tmp_i, current_overlap)
    2606           96 :          DO i = 1, nvec
    2607           76 :             state => wfi_get_snapshot(wf_history, wf_index=i)
    2608           96 :             CALL dbcsr_add(tmp_i, state%overlap, 1.0_dp, -coeffs(i))
    2609              :          END DO
    2610           20 :          error = dbcsr_frobenius_norm(tmp_i)/dbcsr_frobenius_norm(state%overlap)
    2611           20 :          IF (io_unit > 0) THEN
    2612           10 :             WRITE (UNIT=io_unit, FMT="(/,T2,A,F20.10)") "GEXT overlap fitting error:", error
    2613              :          END IF
    2614              :       END IF
    2615              : 
    2616              :       ! free the memory
    2617           20 :       CALL dbcsr_release(tmp_i)
    2618           20 :       CALL dbcsr_release(tmp_j)
    2619           20 :       CALL dbcsr_release(tmp_k)
    2620           20 :       CALL dbcsr_release(target_overlap)
    2621           20 :       DEALLOCATE (A, b)
    2622              : 
    2623           44 :    END SUBROUTINE tr_fitting
    2624              : 
    2625              : END MODULE qs_wf_history_methods
        

Generated by: LCOV version 2.0-1