LCOV - code coverage report
Current view: top level - src - qs_wf_history_methods.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:6d276e9) Lines: 91.1 % 1410 1284
Test Date: 2026-09-10 07:29:18 Functions: 94.7 % 19 18

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

Generated by: LCOV version 2.0-1