LCOV - code coverage report
Current view: top level - src - qs_wf_history_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 425 509 83.5 %
Date: 2024-03-28 07:31:50 Functions: 9 10 90.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 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 cp_control_types,                ONLY: dft_control_type
      27             :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply,&
      28             :                                               dbcsr_allocate_matrix_set,&
      29             :                                               dbcsr_deallocate_matrix_set
      30             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_scale,&
      31             :                                               cp_fm_scale_and_add
      32             :    USE cp_fm_pool_types,                ONLY: cp_fm_pool_p_type,&
      33             :                                               cp_fm_pool_type,&
      34             :                                               fm_pools_create_fm_vect,&
      35             :                                               fm_pools_give_back_fm_vect
      36             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      37             :                                               cp_fm_struct_release,&
      38             :                                               cp_fm_struct_type
      39             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      40             :                                               cp_fm_get_info,&
      41             :                                               cp_fm_release,&
      42             :                                               cp_fm_to_fm,&
      43             :                                               cp_fm_type
      44             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      45             :                                               cp_logger_type,&
      46             :                                               cp_to_string
      47             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      48             :                                               cp_print_key_unit_nr,&
      49             :                                               low_print_level
      50             :    USE dbcsr_api,                       ONLY: dbcsr_add,&
      51             :                                               dbcsr_copy,&
      52             :                                               dbcsr_deallocate_matrix,&
      53             :                                               dbcsr_p_type
      54             :    USE input_constants,                 ONLY: &
      55             :         wfi_aspc_nr, wfi_frozen_method_nr, wfi_linear_p_method_nr, wfi_linear_ps_method_nr, &
      56             :         wfi_linear_wf_method_nr, wfi_ps_method_nr, wfi_use_guess_method_nr, &
      57             :         wfi_use_prev_p_method_nr, wfi_use_prev_rho_r_method_nr, wfi_use_prev_wf_method_nr
      58             :    USE kinds,                           ONLY: dp
      59             :    USE mathlib,                         ONLY: binomial
      60             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      61             :    USE pw_env_types,                    ONLY: pw_env_get,&
      62             :                                               pw_env_type
      63             :    USE pw_methods,                      ONLY: pw_copy
      64             :    USE pw_pool_types,                   ONLY: pw_pool_type
      65             :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      66             :                                               pw_r3d_rs_type
      67             :    USE qs_density_matrices,             ONLY: calculate_density_matrix
      68             :    USE qs_environment_types,            ONLY: get_qs_env,&
      69             :                                               qs_environment_type,&
      70             :                                               set_qs_env
      71             :    USE qs_ks_types,                     ONLY: qs_ks_did_change
      72             :    USE qs_matrix_pools,                 ONLY: mpools_get,&
      73             :                                               qs_matrix_pools_type
      74             :    USE qs_mo_methods,                   ONLY: make_basis_cholesky,&
      75             :                                               make_basis_lowdin,&
      76             :                                               make_basis_simple,&
      77             :                                               make_basis_sm
      78             :    USE qs_mo_types,                     ONLY: get_mo_set,&
      79             :                                               mo_set_type
      80             :    USE qs_rho_methods,                  ONLY: qs_rho_update_rho
      81             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      82             :                                               qs_rho_set,&
      83             :                                               qs_rho_type
      84             :    USE qs_scf_types,                    ONLY: ot_method_nr,&
      85             :                                               qs_scf_env_type
      86             :    USE qs_wf_history_types,             ONLY: qs_wf_history_type,&
      87             :                                               qs_wf_snapshot_type,&
      88             :                                               wfi_get_snapshot,&
      89             :                                               wfi_release
      90             :    USE scf_control_types,               ONLY: scf_control_type
      91             : #include "./base/base_uses.f90"
      92             : 
      93             :    IMPLICIT NONE
      94             :    PRIVATE
      95             : 
      96             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      97             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_wf_history_methods'
      98             : 
      99             :    PUBLIC :: wfi_create, wfi_update, wfi_create_for_kp, &
     100             :              wfi_extrapolate, wfi_get_method_label, &
     101             :              reorthogonalize_vectors, wfi_purge_history
     102             : 
     103             : CONTAINS
     104             : 
     105             : ! **************************************************************************************************
     106             : !> \brief allocates and initialize a wavefunction snapshot
     107             : !> \param snapshot the snapshot to create
     108             : !> \par History
     109             : !>      02.2003 created [fawzi]
     110             : !>      02.2005 added wf_mol [MI]
     111             : !> \author fawzi
     112             : ! **************************************************************************************************
     113        9350 :    SUBROUTINE wfs_create(snapshot)
     114             :       TYPE(qs_wf_snapshot_type), INTENT(OUT)             :: snapshot
     115             : 
     116        9350 :       NULLIFY (snapshot%wf, snapshot%rho_r, &
     117        9350 :                snapshot%rho_g, snapshot%rho_ao, snapshot%rho_ao_kp, &
     118        9350 :                snapshot%overlap, snapshot%rho_frozen)
     119        9350 :       snapshot%dt = 1.0_dp
     120        9350 :    END SUBROUTINE wfs_create
     121             : 
     122             : ! **************************************************************************************************
     123             : !> \brief updates the given snapshot
     124             : !> \param snapshot the snapshot to be updated
     125             : !> \param wf_history the history
     126             : !> \param qs_env the qs_env that should be snapshotted
     127             : !> \param dt the time of the snapshot (wrt. to the previous snapshot)
     128             : !> \par History
     129             : !>      02.2003 created [fawzi]
     130             : !>      02.2005 added kg_fm_mol_set for KG_GPW [MI]
     131             : !> \author fawzi
     132             : ! **************************************************************************************************
     133       16450 :    SUBROUTINE wfs_update(snapshot, wf_history, qs_env, dt)
     134             :       TYPE(qs_wf_snapshot_type), POINTER                 :: snapshot
     135             :       TYPE(qs_wf_history_type), POINTER                  :: wf_history
     136             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     137             :       REAL(KIND=dp), INTENT(in), OPTIONAL                :: dt
     138             : 
     139             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'wfs_update'
     140             : 
     141             :       INTEGER                                            :: handle, img, ispin, nimg, nspins
     142       16450 :       TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER     :: ao_mo_pools
     143             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     144       16450 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s, rho_ao
     145       16450 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: rho_ao_kp
     146             :       TYPE(dft_control_type), POINTER                    :: dft_control
     147       16450 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     148       16450 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g
     149             :       TYPE(pw_env_type), POINTER                         :: pw_env
     150             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     151       16450 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r
     152             :       TYPE(qs_rho_type), POINTER                         :: rho
     153             : 
     154       16450 :       CALL timeset(routineN, handle)
     155             : 
     156       16450 :       NULLIFY (pw_env, auxbas_pw_pool, ao_mo_pools, dft_control, mos, mo_coeff, &
     157       16450 :                rho, rho_r, rho_g, rho_ao, matrix_s)
     158             :       CALL get_qs_env(qs_env, pw_env=pw_env, &
     159       16450 :                       dft_control=dft_control, rho=rho)
     160       16450 :       CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_pools)
     161       16450 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     162             : 
     163       16450 :       CPASSERT(ASSOCIATED(wf_history))
     164       16450 :       CPASSERT(ASSOCIATED(dft_control))
     165       16450 :       IF (.NOT. ASSOCIATED(snapshot)) THEN
     166        9350 :          ALLOCATE (snapshot)
     167        9350 :          CALL wfs_create(snapshot)
     168             :       END IF
     169       16450 :       CPASSERT(wf_history%ref_count > 0)
     170             : 
     171       16450 :       nspins = dft_control%nspins
     172       16450 :       snapshot%dt = 1.0_dp
     173       16450 :       IF (PRESENT(dt)) snapshot%dt = dt
     174       16450 :       IF (wf_history%store_wf) THEN
     175       15876 :          CALL get_qs_env(qs_env, mos=mos)
     176       15876 :          IF (.NOT. ASSOCIATED(snapshot%wf)) THEN
     177             :             CALL fm_pools_create_fm_vect(ao_mo_pools, snapshot%wf, &
     178        9272 :                                          name="ws_snap-ws")
     179        9272 :             CPASSERT(nspins == SIZE(snapshot%wf))
     180             :          END IF
     181       33604 :          DO ispin = 1, nspins
     182       17728 :             CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
     183       33604 :             CALL cp_fm_to_fm(mo_coeff, snapshot%wf(ispin))
     184             :          END DO
     185             :       ELSE
     186         574 :          CALL fm_pools_give_back_fm_vect(ao_mo_pools, snapshot%wf)
     187             :       END IF
     188             : 
     189       16450 :       IF (wf_history%store_rho_r) THEN
     190           0 :          CALL qs_rho_get(rho, rho_r=rho_r)
     191           0 :          CPASSERT(ASSOCIATED(rho_r))
     192           0 :          IF (.NOT. ASSOCIATED(snapshot%rho_r)) THEN
     193           0 :             ALLOCATE (snapshot%rho_r(nspins))
     194           0 :             DO ispin = 1, nspins
     195           0 :                CALL auxbas_pw_pool%create_pw(snapshot%rho_r(ispin))
     196             :             END DO
     197             :          END IF
     198           0 :          DO ispin = 1, nspins
     199           0 :             CALL pw_copy(rho_r(ispin), snapshot%rho_r(ispin))
     200             :          END DO
     201       16450 :       ELSE IF (ASSOCIATED(snapshot%rho_r)) THEN
     202           0 :          DO ispin = 1, SIZE(snapshot%rho_r)
     203           0 :             CALL auxbas_pw_pool%give_back_pw(snapshot%rho_r(ispin))
     204             :          END DO
     205           0 :          DEALLOCATE (snapshot%rho_r)
     206             :       END IF
     207             : 
     208       16450 :       IF (wf_history%store_rho_g) THEN
     209           0 :          CALL qs_rho_get(rho, rho_g=rho_g)
     210           0 :          CPASSERT(ASSOCIATED(rho_g))
     211           0 :          IF (.NOT. ASSOCIATED(snapshot%rho_g)) THEN
     212           0 :             ALLOCATE (snapshot%rho_g(nspins))
     213           0 :             DO ispin = 1, nspins
     214           0 :                CALL auxbas_pw_pool%create_pw(snapshot%rho_g(ispin))
     215             :             END DO
     216             :          END IF
     217           0 :          DO ispin = 1, nspins
     218           0 :             CALL pw_copy(rho_g(ispin), snapshot%rho_g(ispin))
     219             :          END DO
     220       16450 :       ELSE IF (ASSOCIATED(snapshot%rho_g)) THEN
     221           0 :          DO ispin = 1, SIZE(snapshot%rho_g)
     222           0 :             CALL auxbas_pw_pool%give_back_pw(snapshot%rho_g(ispin))
     223             :          END DO
     224           0 :          DEALLOCATE (snapshot%rho_g)
     225             :       END IF
     226             : 
     227       16450 :       IF (ASSOCIATED(snapshot%rho_ao)) THEN ! the sparsity might be different
     228             :          ! (future struct:check)
     229         270 :          CALL dbcsr_deallocate_matrix_set(snapshot%rho_ao)
     230             :       END IF
     231       16450 :       IF (wf_history%store_rho_ao) THEN
     232         338 :          CALL qs_rho_get(rho, rho_ao=rho_ao)
     233         338 :          CPASSERT(ASSOCIATED(rho_ao))
     234             : 
     235         338 :          CALL dbcsr_allocate_matrix_set(snapshot%rho_ao, nspins)
     236         836 :          DO ispin = 1, nspins
     237         498 :             ALLOCATE (snapshot%rho_ao(ispin)%matrix)
     238         836 :             CALL dbcsr_copy(snapshot%rho_ao(ispin)%matrix, rho_ao(ispin)%matrix)
     239             :          END DO
     240             :       END IF
     241             : 
     242       16450 :       IF (ASSOCIATED(snapshot%rho_ao_kp)) THEN ! the sparsity might be different
     243             :          ! (future struct:check)
     244         214 :          CALL dbcsr_deallocate_matrix_set(snapshot%rho_ao_kp)
     245             :       END IF
     246       16450 :       IF (wf_history%store_rho_ao_kp) THEN
     247         220 :          CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
     248         220 :          CPASSERT(ASSOCIATED(rho_ao_kp))
     249             : 
     250         220 :          nimg = dft_control%nimages
     251         220 :          CALL dbcsr_allocate_matrix_set(snapshot%rho_ao_kp, nspins, nimg)
     252         530 :          DO ispin = 1, nspins
     253       31010 :             DO img = 1, nimg
     254       30480 :                ALLOCATE (snapshot%rho_ao_kp(ispin, img)%matrix)
     255             :                CALL dbcsr_copy(snapshot%rho_ao_kp(ispin, img)%matrix, &
     256       30790 :                                rho_ao_kp(ispin, img)%matrix)
     257             :             END DO
     258             :          END DO
     259             :       END IF
     260             : 
     261       16450 :       IF (ASSOCIATED(snapshot%overlap)) THEN ! the sparsity might be different
     262             :          ! (future struct:check)
     263        3936 :          CALL dbcsr_deallocate_matrix(snapshot%overlap)
     264             :       END IF
     265       16450 :       IF (wf_history%store_overlap) THEN
     266       12184 :          CALL get_qs_env(qs_env, matrix_s=matrix_s)
     267       12184 :          CPASSERT(ASSOCIATED(matrix_s))
     268       12184 :          CPASSERT(ASSOCIATED(matrix_s(1)%matrix))
     269       12184 :          ALLOCATE (snapshot%overlap)
     270       12184 :          CALL dbcsr_copy(snapshot%overlap, matrix_s(1)%matrix)
     271             :       END IF
     272             : 
     273             :       IF (wf_history%store_frozen_density) THEN
     274             :          ! do nothing
     275             :          ! CALL deallocate_matrix_set(snapshot%rho_frozen%rho_ao)
     276             :       END IF
     277             : 
     278       16450 :       CALL timestop(handle)
     279             : 
     280       16450 :    END SUBROUTINE wfs_update
     281             : 
     282             : ! **************************************************************************************************
     283             : !> \brief ...
     284             : !> \param wf_history ...
     285             : !> \param interpolation_method_nr the tag of the method used for
     286             : !>        the extrapolation of the initial density for the next md step
     287             : !>        (see qs_wf_history_types:wfi_*_method_nr)
     288             : !> \param extrapolation_order ...
     289             : !> \param has_unit_metric ...
     290             : !> \par History
     291             : !>      02.2003 created [fawzi]
     292             : !> \author fawzi
     293             : ! **************************************************************************************************
     294        6556 :    SUBROUTINE wfi_create(wf_history, interpolation_method_nr, extrapolation_order, &
     295             :                          has_unit_metric)
     296             :       TYPE(qs_wf_history_type), POINTER                  :: wf_history
     297             :       INTEGER, INTENT(in)                                :: interpolation_method_nr, &
     298             :                                                             extrapolation_order
     299             :       LOGICAL, INTENT(IN)                                :: has_unit_metric
     300             : 
     301             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'wfi_create'
     302             : 
     303             :       INTEGER                                            :: handle, i
     304             : 
     305        6556 :       CALL timeset(routineN, handle)
     306             : 
     307        6556 :       ALLOCATE (wf_history)
     308        6556 :       wf_history%ref_count = 1
     309        6556 :       wf_history%memory_depth = 0
     310        6556 :       wf_history%snapshot_count = 0
     311        6556 :       wf_history%last_state_index = 1
     312        6556 :       wf_history%store_wf = .FALSE.
     313        6556 :       wf_history%store_rho_r = .FALSE.
     314        6556 :       wf_history%store_rho_g = .FALSE.
     315        6556 :       wf_history%store_rho_ao = .FALSE.
     316        6556 :       wf_history%store_rho_ao_kp = .FALSE.
     317        6556 :       wf_history%store_overlap = .FALSE.
     318        6556 :       wf_history%store_frozen_density = .FALSE.
     319        6556 :       NULLIFY (wf_history%past_states)
     320             : 
     321        6556 :       wf_history%interpolation_method_nr = interpolation_method_nr
     322             : 
     323             :       SELECT CASE (wf_history%interpolation_method_nr)
     324             :       CASE (wfi_use_guess_method_nr)
     325             :          wf_history%memory_depth = 0
     326             :       CASE (wfi_use_prev_wf_method_nr)
     327          62 :          wf_history%memory_depth = 0
     328             :       CASE (wfi_use_prev_p_method_nr)
     329          62 :          wf_history%memory_depth = 1
     330          62 :          wf_history%store_rho_ao = .TRUE.
     331             :       CASE (wfi_use_prev_rho_r_method_nr)
     332           4 :          wf_history%memory_depth = 1
     333           4 :          wf_history%store_rho_ao = .TRUE.
     334             :       CASE (wfi_linear_wf_method_nr)
     335           4 :          wf_history%memory_depth = 2
     336           4 :          wf_history%store_wf = .TRUE.
     337             :       CASE (wfi_linear_p_method_nr)
     338           4 :          wf_history%memory_depth = 2
     339           4 :          wf_history%store_rho_ao = .TRUE.
     340             :       CASE (wfi_linear_ps_method_nr)
     341           6 :          wf_history%memory_depth = 2
     342           6 :          wf_history%store_wf = .TRUE.
     343           6 :          IF (.NOT. has_unit_metric) wf_history%store_overlap = .TRUE.
     344             :       CASE (wfi_ps_method_nr)
     345         337 :          CALL cite_reference(VandeVondele2005a)
     346         337 :          wf_history%memory_depth = extrapolation_order + 1
     347         337 :          wf_history%store_wf = .TRUE.
     348         337 :          IF (.NOT. has_unit_metric) wf_history%store_overlap = .TRUE.
     349             :       CASE (wfi_frozen_method_nr)
     350           4 :          wf_history%memory_depth = 1
     351           4 :          wf_history%store_frozen_density = .TRUE.
     352             :       CASE (wfi_aspc_nr)
     353        5939 :          wf_history%memory_depth = extrapolation_order + 2
     354        5939 :          wf_history%store_wf = .TRUE.
     355        5939 :          IF (.NOT. has_unit_metric) wf_history%store_overlap = .TRUE.
     356             :       CASE default
     357             :          CALL cp_abort(__LOCATION__, &
     358             :                        "Unknown interpolation method: "// &
     359           0 :                        TRIM(ADJUSTL(cp_to_string(interpolation_method_nr))))
     360        6556 :          wf_history%interpolation_method_nr = wfi_use_prev_rho_r_method_nr
     361             :       END SELECT
     362       19472 :       ALLOCATE (wf_history%past_states(wf_history%memory_depth))
     363             : 
     364       37379 :       DO i = 1, SIZE(wf_history%past_states)
     365       37379 :          NULLIFY (wf_history%past_states(i)%snapshot)
     366             :       END DO
     367             : 
     368        6556 :       CALL timestop(handle)
     369        6556 :    END SUBROUTINE wfi_create
     370             : 
     371             : ! **************************************************************************************************
     372             : !> \brief ...
     373             : !> \param wf_history ...
     374             : !> \par History
     375             : !>      06.2015 created [jhu]
     376             : !> \author fawzi
     377             : ! **************************************************************************************************
     378         148 :    SUBROUTINE wfi_create_for_kp(wf_history)
     379             :       TYPE(qs_wf_history_type), POINTER                  :: wf_history
     380             : 
     381         148 :       CPASSERT(ASSOCIATED(wf_history))
     382         148 :       IF (wf_history%store_rho_ao) THEN
     383           6 :          wf_history%store_rho_ao_kp = .TRUE.
     384           6 :          wf_history%store_rho_ao = .FALSE.
     385             :       END IF
     386             :       ! Check for KP compatible methods
     387         148 :       IF (wf_history%store_wf) THEN
     388           0 :          CPABORT("WFN based interpolation method not possible for kpoints.")
     389             :       END IF
     390         148 :       IF (wf_history%store_frozen_density) THEN
     391           0 :          CPABORT("Frozen density initialization method not possible for kpoints.")
     392             :       END IF
     393         148 :       IF (wf_history%store_overlap) THEN
     394           0 :          CPABORT("Inconsistent interpolation method for kpoints.")
     395             :       END IF
     396             : 
     397         148 :    END SUBROUTINE wfi_create_for_kp
     398             : 
     399             : ! **************************************************************************************************
     400             : !> \brief returns a string describing the interpolation method
     401             : !> \param method_nr ...
     402             : !> \return ...
     403             : !> \par History
     404             : !>      02.2003 created [fawzi]
     405             : !> \author fawzi
     406             : ! **************************************************************************************************
     407        8930 :    FUNCTION wfi_get_method_label(method_nr) RESULT(res)
     408             :       INTEGER, INTENT(in)                                :: method_nr
     409             :       CHARACTER(len=30)                                  :: res
     410             : 
     411        8930 :       res = "unknown"
     412        9166 :       SELECT CASE (method_nr)
     413             :       CASE (wfi_use_prev_p_method_nr)
     414         236 :          res = "previous_p"
     415             :       CASE (wfi_use_prev_wf_method_nr)
     416         210 :          res = "previous_wf"
     417             :       CASE (wfi_use_prev_rho_r_method_nr)
     418           4 :          res = "previous_rho_r"
     419             :       CASE (wfi_use_guess_method_nr)
     420        2960 :          res = "initial_guess"
     421             :       CASE (wfi_linear_wf_method_nr)
     422           2 :          res = "mo linear"
     423             :       CASE (wfi_linear_p_method_nr)
     424           2 :          res = "P linear"
     425             :       CASE (wfi_linear_ps_method_nr)
     426           6 :          res = "PS linear"
     427             :       CASE (wfi_ps_method_nr)
     428         124 :          res = "PS Nth order"
     429             :       CASE (wfi_frozen_method_nr)
     430           4 :          res = "frozen density approximation"
     431             :       CASE (wfi_aspc_nr)
     432        5382 :          res = "ASPC"
     433             :       CASE default
     434             :          CALL cp_abort(__LOCATION__, &
     435             :                        "Unknown interpolation method: "// &
     436        8930 :                        TRIM(ADJUSTL(cp_to_string(method_nr))))
     437             :       END SELECT
     438        8930 :    END FUNCTION wfi_get_method_label
     439             : 
     440             : ! **************************************************************************************************
     441             : !> \brief calculates the new starting state for the scf for the next
     442             : !>      wf optimization
     443             : !> \param wf_history the previous history needed to extrapolate
     444             : !> \param qs_env the qs env with the latest result, and that will contain
     445             : !>        the new starting state
     446             : !> \param dt the time at which to extrapolate (wrt. to the last snapshot)
     447             : !> \param extrapolation_method_nr returns the extrapolation method used
     448             : !> \param orthogonal_wf ...
     449             : !> \par History
     450             : !>      02.2003 created [fawzi]
     451             : !>      11.2003 Joost VandeVondele : Implemented Nth order PS extrapolation
     452             : !> \author fawzi
     453             : ! **************************************************************************************************
     454       17495 :    SUBROUTINE wfi_extrapolate(wf_history, qs_env, dt, extrapolation_method_nr, &
     455             :                               orthogonal_wf)
     456             :       TYPE(qs_wf_history_type), POINTER                  :: wf_history
     457             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     458             :       REAL(KIND=dp), INTENT(IN)                          :: dt
     459             :       INTEGER, INTENT(OUT), OPTIONAL                     :: extrapolation_method_nr
     460             :       LOGICAL, INTENT(OUT), OPTIONAL                     :: orthogonal_wf
     461             : 
     462             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'wfi_extrapolate'
     463             : 
     464             :       INTEGER                                            :: actual_extrapolation_method_nr, handle, &
     465             :                                                             i, img, io_unit, ispin, k, n, nmo, &
     466             :                                                             nvec, print_level
     467             :       LOGICAL                                            :: do_kpoints, my_orthogonal_wf, use_overlap
     468             :       REAL(KIND=dp)                                      :: alpha, t0, t1, t2
     469       17495 :       TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER     :: ao_mo_fm_pools
     470             :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct, matrix_struct_new
     471             :       TYPE(cp_fm_type)                                   :: csc, fm_tmp
     472             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     473             :       TYPE(cp_logger_type), POINTER                      :: logger
     474       17495 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao, rho_frozen_ao
     475       17495 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: rho_ao_kp
     476       17495 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     477             :       TYPE(qs_rho_type), POINTER                         :: rho
     478             :       TYPE(qs_wf_snapshot_type), POINTER                 :: t0_state, t1_state
     479             : 
     480       17495 :       NULLIFY (mos, ao_mo_fm_pools, t0_state, t1_state, mo_coeff, &
     481       17495 :                rho, rho_ao, rho_frozen_ao)
     482             : 
     483       17495 :       use_overlap = wf_history%store_overlap
     484             : 
     485       17495 :       CALL timeset(routineN, handle)
     486       17495 :       logger => cp_get_default_logger()
     487       17495 :       print_level = logger%iter_info%print_level
     488             :       io_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%SCF%PRINT%PROGRAM_RUN_INFO", &
     489       17495 :                                      extension=".scfLog")
     490             : 
     491       17495 :       CPASSERT(ASSOCIATED(wf_history))
     492       17495 :       CPASSERT(wf_history%ref_count > 0)
     493       17495 :       CPASSERT(ASSOCIATED(qs_env))
     494       17495 :       CALL get_qs_env(qs_env, mos=mos, rho=rho, do_kpoints=do_kpoints)
     495       17495 :       CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
     496             :       ! chooses the method for this extrapolation
     497       17495 :       IF (wf_history%snapshot_count < 1) THEN
     498             :          actual_extrapolation_method_nr = wfi_use_guess_method_nr
     499             :       ELSE
     500       12088 :          actual_extrapolation_method_nr = wf_history%interpolation_method_nr
     501             :       END IF
     502             : 
     503           8 :       SELECT CASE (actual_extrapolation_method_nr)
     504             :       CASE (wfi_linear_wf_method_nr)
     505           8 :          IF (wf_history%snapshot_count < 2) THEN
     506           4 :             actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
     507             :          END IF
     508             :       CASE (wfi_linear_p_method_nr)
     509           8 :          IF (wf_history%snapshot_count < 2) THEN
     510           4 :             IF (do_kpoints) THEN
     511             :                actual_extrapolation_method_nr = wfi_use_guess_method_nr
     512             :             ELSE
     513           4 :                actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
     514             :             END IF
     515             :          END IF
     516             :       CASE (wfi_linear_ps_method_nr)
     517       12088 :          IF (wf_history%snapshot_count < 2) THEN
     518           6 :             actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
     519             :          END IF
     520             :       END SELECT
     521             : 
     522       17495 :       IF (PRESENT(extrapolation_method_nr)) &
     523       17495 :          extrapolation_method_nr = actual_extrapolation_method_nr
     524       17495 :       my_orthogonal_wf = .FALSE.
     525             : 
     526           8 :       SELECT CASE (actual_extrapolation_method_nr)
     527             :       CASE (wfi_frozen_method_nr)
     528           8 :          CPASSERT(.NOT. do_kpoints)
     529           8 :          t0_state => wfi_get_snapshot(wf_history, wf_index=1)
     530           8 :          CPASSERT(ASSOCIATED(t0_state%rho_frozen))
     531             : 
     532           8 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     533           8 :          CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
     534             : 
     535           8 :          CALL qs_rho_get(t0_state%rho_frozen, rho_ao=rho_frozen_ao)
     536           8 :          CALL qs_rho_get(rho, rho_ao=rho_ao)
     537          16 :          DO ispin = 1, SIZE(rho_frozen_ao)
     538             :             CALL dbcsr_copy(rho_ao(ispin)%matrix, &
     539             :                             rho_frozen_ao(ispin)%matrix, &
     540          16 :                             keep_sparsity=.TRUE.)
     541             :          END DO
     542             :          !FM updating rho_ao directly with t0_state%rho_ao would have the
     543             :          !FM wrong matrix structure
     544           8 :          CALL qs_rho_update_rho(rho, qs_env=qs_env)
     545           8 :          CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
     546             : 
     547           8 :          my_orthogonal_wf = .FALSE.
     548             :       CASE (wfi_use_prev_rho_r_method_nr)
     549           8 :          t0_state => wfi_get_snapshot(wf_history, wf_index=1)
     550           8 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     551           8 :          CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
     552           8 :          IF (do_kpoints) THEN
     553           0 :             CPASSERT(ASSOCIATED(t0_state%rho_ao_kp))
     554           0 :             CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
     555           0 :             DO ispin = 1, SIZE(t0_state%rho_ao_kp, 1)
     556           0 :                DO img = 1, SIZE(t0_state%rho_ao_kp, 2)
     557           0 :                   IF (img > SIZE(rho_ao_kp, 2)) THEN
     558           0 :                      CPWARN("Change in cell neighborlist: might affect quality of initial guess")
     559             :                   ELSE
     560             :                      CALL dbcsr_copy(rho_ao_kp(ispin, img)%matrix, &
     561             :                                      t0_state%rho_ao_kp(ispin, img)%matrix, &
     562           0 :                                      keep_sparsity=.TRUE.)
     563             :                   END IF
     564             :                END DO
     565             :             END DO
     566             :          ELSE
     567           8 :             CPASSERT(ASSOCIATED(t0_state%rho_ao))
     568           8 :             CALL qs_rho_get(rho, rho_ao=rho_ao)
     569          16 :             DO ispin = 1, SIZE(t0_state%rho_ao)
     570             :                CALL dbcsr_copy(rho_ao(ispin)%matrix, &
     571             :                                t0_state%rho_ao(ispin)%matrix, &
     572          16 :                                keep_sparsity=.TRUE.)
     573             :             END DO
     574             :          END IF
     575             :          ! Why is rho_g valid at this point ?
     576           8 :          CALL qs_rho_set(rho, rho_g_valid=.TRUE.)
     577             : 
     578             :          ! does nothing
     579             :       CASE (wfi_use_prev_wf_method_nr)
     580         420 :          CPASSERT(.NOT. do_kpoints)
     581         420 :          my_orthogonal_wf = .TRUE.
     582         420 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     583         420 :          CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
     584         420 :          CALL qs_rho_get(rho, rho_ao=rho_ao)
     585         888 :          DO ispin = 1, SIZE(mos)
     586             :             CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, &
     587         468 :                             nmo=nmo)
     588             :             CALL reorthogonalize_vectors(qs_env, &
     589             :                                          v_matrix=mo_coeff, &
     590         468 :                                          n_col=nmo)
     591             :             CALL calculate_density_matrix(mo_set=mos(ispin), &
     592        1356 :                                           density_matrix=rho_ao(ispin)%matrix)
     593             :          END DO
     594         420 :          CALL qs_rho_update_rho(rho, qs_env=qs_env)
     595         420 :          CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
     596             : 
     597             :       CASE (wfi_use_prev_p_method_nr)
     598         472 :          t0_state => wfi_get_snapshot(wf_history, wf_index=1)
     599         472 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     600         472 :          CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
     601         472 :          IF (do_kpoints) THEN
     602         214 :             CPASSERT(ASSOCIATED(t0_state%rho_ao_kp))
     603         214 :             CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
     604         516 :             DO ispin = 1, SIZE(t0_state%rho_ao_kp, 1)
     605       30222 :                DO img = 1, SIZE(t0_state%rho_ao_kp, 2)
     606       30008 :                   IF (img > SIZE(rho_ao_kp, 2)) THEN
     607           0 :                      CPWARN("Change in cell neighborlist: might affect quality of initial guess")
     608             :                   ELSE
     609             :                      CALL dbcsr_copy(rho_ao_kp(ispin, img)%matrix, &
     610             :                                      t0_state%rho_ao_kp(ispin, img)%matrix, &
     611       29706 :                                      keep_sparsity=.TRUE.)
     612             :                   END IF
     613             :                END DO
     614             :             END DO
     615             :          ELSE
     616         258 :             CPASSERT(ASSOCIATED(t0_state%rho_ao))
     617         258 :             CALL qs_rho_get(rho, rho_ao=rho_ao)
     618         646 :             DO ispin = 1, SIZE(t0_state%rho_ao)
     619             :                CALL dbcsr_copy(rho_ao(ispin)%matrix, &
     620             :                                t0_state%rho_ao(ispin)%matrix, &
     621         646 :                                keep_sparsity=.TRUE.)
     622             :             END DO
     623             :          END IF
     624             :          !FM updating rho_ao directly with t0_state%rho_ao would have the
     625             :          !FM wrong matrix structure
     626         472 :          CALL qs_rho_update_rho(rho, qs_env=qs_env)
     627         472 :          CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
     628             : 
     629             :       CASE (wfi_use_guess_method_nr)
     630             :          !FM more clean to do it here, but it
     631             :          !FM might need to read a file (restart) and thus globenv
     632             :          !FM I do not want globenv here, thus done by the caller
     633             :          !FM (btw. it also needs the eigensolver, and unless you relocate it
     634             :          !FM gives circular dependencies)
     635        5859 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     636        5859 :          CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
     637             :       CASE (wfi_linear_wf_method_nr)
     638           4 :          CPASSERT(.NOT. do_kpoints)
     639           4 :          t0_state => wfi_get_snapshot(wf_history, wf_index=2)
     640           4 :          t1_state => wfi_get_snapshot(wf_history, wf_index=1)
     641           4 :          CPASSERT(ASSOCIATED(t0_state))
     642           4 :          CPASSERT(ASSOCIATED(t1_state))
     643           4 :          CPASSERT(ASSOCIATED(t0_state%wf))
     644           4 :          CPASSERT(ASSOCIATED(t1_state%wf))
     645           4 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     646           4 :          CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
     647             : 
     648           4 :          my_orthogonal_wf = .TRUE.
     649           4 :          t0 = 0.0_dp
     650           4 :          t1 = t1_state%dt
     651           4 :          t2 = t1 + dt
     652           4 :          CALL qs_rho_get(rho, rho_ao=rho_ao)
     653           8 :          DO ispin = 1, SIZE(mos)
     654             :             CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, &
     655           4 :                             nmo=nmo)
     656             :             CALL cp_fm_scale_and_add(alpha=0.0_dp, &
     657             :                                      matrix_a=mo_coeff, &
     658             :                                      matrix_b=t1_state%wf(ispin), &
     659           4 :                                      beta=(t2 - t0)/(t1 - t0))
     660             :             ! this copy should be unnecessary
     661             :             CALL cp_fm_scale_and_add(alpha=1.0_dp, &
     662             :                                      matrix_a=mo_coeff, &
     663           4 :                                      beta=(t1 - t2)/(t1 - t0), matrix_b=t0_state%wf(ispin))
     664             :             CALL reorthogonalize_vectors(qs_env, &
     665             :                                          v_matrix=mo_coeff, &
     666           4 :                                          n_col=nmo)
     667             :             CALL calculate_density_matrix(mo_set=mos(ispin), &
     668          12 :                                           density_matrix=rho_ao(ispin)%matrix)
     669             :          END DO
     670           4 :          CALL qs_rho_update_rho(rho, qs_env=qs_env)
     671             : 
     672             :          CALL qs_ks_did_change(qs_env%ks_env, &
     673           4 :                                rho_changed=.TRUE.)
     674             :       CASE (wfi_linear_p_method_nr)
     675           4 :          t0_state => wfi_get_snapshot(wf_history, wf_index=2)
     676           4 :          t1_state => wfi_get_snapshot(wf_history, wf_index=1)
     677           4 :          CPASSERT(ASSOCIATED(t0_state))
     678           4 :          CPASSERT(ASSOCIATED(t1_state))
     679           4 :          IF (do_kpoints) THEN
     680           0 :             CPASSERT(ASSOCIATED(t0_state%rho_ao_kp))
     681           0 :             CPASSERT(ASSOCIATED(t1_state%rho_ao_kp))
     682             :          ELSE
     683           4 :             CPASSERT(ASSOCIATED(t0_state%rho_ao))
     684           4 :             CPASSERT(ASSOCIATED(t1_state%rho_ao))
     685             :          END IF
     686           4 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     687           4 :          CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
     688             : 
     689           4 :          t0 = 0.0_dp
     690           4 :          t1 = t1_state%dt
     691           4 :          t2 = t1 + dt
     692           4 :          IF (do_kpoints) THEN
     693           0 :             CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
     694           0 :             DO ispin = 1, SIZE(rho_ao_kp, 1)
     695           0 :                DO img = 1, SIZE(rho_ao_kp, 2)
     696           0 :                   IF (img > SIZE(t0_state%rho_ao_kp, 2) .OR. &
     697           0 :                       img > SIZE(t1_state%rho_ao_kp, 2)) THEN
     698           0 :                      CPWARN("Change in cell neighborlist: might affect quality of initial guess")
     699             :                   ELSE
     700             :                      CALL dbcsr_add(rho_ao_kp(ispin, img)%matrix, t1_state%rho_ao_kp(ispin, img)%matrix, &
     701           0 :                                     alpha_scalar=0.0_dp, beta_scalar=(t2 - t0)/(t1 - t0)) ! this copy should be unnecessary
     702             :                      CALL dbcsr_add(rho_ao_kp(ispin, img)%matrix, t0_state%rho_ao_kp(ispin, img)%matrix, &
     703           0 :                                     alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0))
     704             :                   END IF
     705             :                END DO
     706             :             END DO
     707             :          ELSE
     708           4 :             CALL qs_rho_get(rho, rho_ao=rho_ao)
     709           8 :             DO ispin = 1, SIZE(rho_ao)
     710             :                CALL dbcsr_add(rho_ao(ispin)%matrix, t1_state%rho_ao(ispin)%matrix, &
     711           4 :                               alpha_scalar=0.0_dp, beta_scalar=(t2 - t0)/(t1 - t0)) ! this copy should be unnecessary
     712             :                CALL dbcsr_add(rho_ao(ispin)%matrix, t0_state%rho_ao(ispin)%matrix, &
     713           8 :                               alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0))
     714             :             END DO
     715             :          END IF
     716           4 :          CALL qs_rho_update_rho(rho, qs_env=qs_env)
     717           4 :          CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
     718             :       CASE (wfi_linear_ps_method_nr)
     719             :          ! wf not calculated, extract with PSC renormalized?
     720             :          ! use wf_linear?
     721          12 :          CPASSERT(.NOT. do_kpoints)
     722          12 :          t0_state => wfi_get_snapshot(wf_history, wf_index=2)
     723          12 :          t1_state => wfi_get_snapshot(wf_history, wf_index=1)
     724          12 :          CPASSERT(ASSOCIATED(t0_state))
     725          12 :          CPASSERT(ASSOCIATED(t1_state))
     726          12 :          CPASSERT(ASSOCIATED(t0_state%wf))
     727          12 :          CPASSERT(ASSOCIATED(t1_state%wf))
     728          12 :          IF (wf_history%store_overlap) THEN
     729           4 :             CPASSERT(ASSOCIATED(t0_state%overlap))
     730           4 :             CPASSERT(ASSOCIATED(t1_state%overlap))
     731             :          END IF
     732          12 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     733          12 :          IF (nvec >= wf_history%memory_depth) THEN
     734          12 :             IF ((qs_env%scf_control%max_scf_hist .NE. 0) .AND. (qs_env%scf_control%eps_scf_hist .NE. 0)) THEN
     735           0 :                qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
     736           0 :                qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
     737           0 :                qs_env%scf_control%outer_scf%have_scf = .FALSE.
     738          12 :             ELSE IF (qs_env%scf_control%max_scf_hist .NE. 0) THEN
     739           0 :                qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
     740           0 :                qs_env%scf_control%outer_scf%have_scf = .FALSE.
     741          12 :             ELSE IF (qs_env%scf_control%eps_scf_hist .NE. 0) THEN
     742           0 :                qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
     743             :             END IF
     744             :          END IF
     745             : 
     746          12 :          my_orthogonal_wf = .TRUE.
     747             :          ! use PS_2=2 PS_1-PS_0
     748             :          ! C_2 comes from using PS_2 as a projector acting on C_1
     749          12 :          CALL qs_rho_get(rho, rho_ao=rho_ao)
     750          24 :          DO ispin = 1, SIZE(mos)
     751          12 :             NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
     752          12 :             CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
     753             :             CALL cp_fm_get_info(mo_coeff, nrow_global=n, ncol_global=k, &
     754          12 :                                 matrix_struct=matrix_struct)
     755             :             CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, &
     756          12 :                                      nrow_global=k, ncol_global=k)
     757          12 :             CALL cp_fm_create(csc, matrix_struct_new)
     758          12 :             CALL cp_fm_struct_release(matrix_struct_new)
     759             : 
     760          12 :             IF (use_overlap) THEN
     761           4 :                CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), mo_coeff, k)
     762           4 :                CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), mo_coeff, 0.0_dp, csc)
     763             :             ELSE
     764             :                CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
     765           8 :                                   t1_state%wf(ispin), 0.0_dp, csc)
     766             :             END IF
     767          12 :             CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, mo_coeff)
     768          12 :             CALL cp_fm_release(csc)
     769          12 :             CALL cp_fm_scale_and_add(-1.0_dp, mo_coeff, 2.0_dp, t1_state%wf(ispin))
     770             :             CALL reorthogonalize_vectors(qs_env, &
     771             :                                          v_matrix=mo_coeff, &
     772          12 :                                          n_col=k)
     773             :             CALL calculate_density_matrix(mo_set=mos(ispin), &
     774          36 :                                           density_matrix=rho_ao(ispin)%matrix)
     775             :          END DO
     776          12 :          CALL qs_rho_update_rho(rho, qs_env=qs_env)
     777          12 :          CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
     778             : 
     779             :       CASE (wfi_ps_method_nr)
     780         248 :          CPASSERT(.NOT. do_kpoints)
     781             :          ! figure out the actual number of vectors to use in the extrapolation:
     782         248 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     783         248 :          CPASSERT(nvec .GT. 0)
     784         248 :          IF (nvec >= wf_history%memory_depth) THEN
     785          58 :             IF ((qs_env%scf_control%max_scf_hist .NE. 0) .AND. (qs_env%scf_control%eps_scf_hist .NE. 0)) THEN
     786           0 :                qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
     787           0 :                qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
     788           0 :                qs_env%scf_control%outer_scf%have_scf = .FALSE.
     789          58 :             ELSE IF (qs_env%scf_control%max_scf_hist .NE. 0) THEN
     790           0 :                qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
     791           0 :                qs_env%scf_control%outer_scf%have_scf = .FALSE.
     792          58 :             ELSE IF (qs_env%scf_control%eps_scf_hist .NE. 0) THEN
     793           0 :                qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
     794             :             END IF
     795             :          END IF
     796             : 
     797         248 :          my_orthogonal_wf = .TRUE.
     798         574 :          DO ispin = 1, SIZE(mos)
     799         326 :             NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
     800         326 :             CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
     801             :             CALL cp_fm_get_info(mo_coeff, nrow_global=n, ncol_global=k, &
     802         326 :                                 matrix_struct=matrix_struct)
     803         326 :             CALL cp_fm_create(fm_tmp, matrix_struct)
     804             :             CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, &
     805         326 :                                      nrow_global=k, ncol_global=k)
     806         326 :             CALL cp_fm_create(csc, matrix_struct_new)
     807         326 :             CALL cp_fm_struct_release(matrix_struct_new)
     808             :             ! first the most recent
     809         326 :             t1_state => wfi_get_snapshot(wf_history, wf_index=1)
     810         326 :             CALL cp_fm_to_fm(t1_state%wf(ispin), mo_coeff)
     811         326 :             alpha = nvec
     812         326 :             CALL cp_fm_scale(alpha, mo_coeff)
     813         326 :             CALL qs_rho_get(rho, rho_ao=rho_ao)
     814         596 :             DO i = 2, nvec
     815         270 :                t0_state => wfi_get_snapshot(wf_history, wf_index=i)
     816         270 :                IF (use_overlap) THEN
     817         232 :                   CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
     818         232 :                   CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
     819             :                ELSE
     820             :                   CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
     821          38 :                                      t1_state%wf(ispin), 0.0_dp, csc)
     822             :                END IF
     823         270 :                CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
     824         270 :                alpha = -1.0_dp*alpha*REAL(nvec - i + 1, dp)/REAL(i, dp)
     825         596 :                CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp)
     826             :             END DO
     827             : 
     828         326 :             CALL cp_fm_release(csc)
     829         326 :             CALL cp_fm_release(fm_tmp)
     830             :             CALL reorthogonalize_vectors(qs_env, &
     831             :                                          v_matrix=mo_coeff, &
     832         326 :                                          n_col=k)
     833             :             CALL calculate_density_matrix(mo_set=mos(ispin), &
     834         900 :                                           density_matrix=rho_ao(ispin)%matrix)
     835             :          END DO
     836         248 :          CALL qs_rho_update_rho(rho, qs_env=qs_env)
     837         248 :          CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
     838             : 
     839             :       CASE (wfi_aspc_nr)
     840       10460 :          CPASSERT(.NOT. do_kpoints)
     841       10460 :          CALL cite_reference(Kolafa2004)
     842       10460 :          CALL cite_reference(Kuhne2007)
     843             :          ! figure out the actual number of vectors to use in the extrapolation:
     844       10460 :          nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
     845       10460 :          CPASSERT(nvec .GT. 0)
     846       10460 :          IF (nvec >= wf_history%memory_depth) THEN
     847        6530 :             IF ((qs_env%scf_control%max_scf_hist .NE. 0) .AND. &
     848             :                 (qs_env%scf_control%eps_scf_hist .NE. 0)) THEN
     849          18 :                qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
     850          18 :                qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
     851          18 :                qs_env%scf_control%outer_scf%have_scf = .FALSE.
     852        6512 :             ELSE IF (qs_env%scf_control%max_scf_hist .NE. 0) THEN
     853          62 :                qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
     854          62 :                qs_env%scf_control%outer_scf%have_scf = .FALSE.
     855        6450 :             ELSE IF (qs_env%scf_control%eps_scf_hist .NE. 0) THEN
     856           8 :                qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
     857             :             END IF
     858             :          END IF
     859             : 
     860       10460 :          my_orthogonal_wf = .TRUE.
     861       10460 :          CALL qs_rho_get(rho, rho_ao=rho_ao)
     862       21381 :          DO ispin = 1, SIZE(mos)
     863       10921 :             NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
     864       10921 :             CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
     865             :             CALL cp_fm_get_info(mo_coeff, &
     866             :                                 nrow_global=n, &
     867             :                                 ncol_global=k, &
     868       10921 :                                 matrix_struct=matrix_struct)
     869       10921 :             CALL cp_fm_create(fm_tmp, matrix_struct)
     870             :             CALL cp_fm_struct_create(matrix_struct_new, &
     871             :                                      template_fmstruct=matrix_struct, &
     872             :                                      nrow_global=k, &
     873       10921 :                                      ncol_global=k)
     874       10921 :             CALL cp_fm_create(csc, matrix_struct_new)
     875       10921 :             CALL cp_fm_struct_release(matrix_struct_new)
     876             :             ! first the most recent
     877             :             t1_state => wfi_get_snapshot(wf_history, &
     878       10921 :                                          wf_index=1)
     879       10921 :             CALL cp_fm_to_fm(t1_state%wf(ispin), mo_coeff)
     880       10921 :             alpha = REAL(4*nvec - 2, KIND=dp)/REAL(nvec + 1, KIND=dp)
     881       10921 :             IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
     882             :                WRITE (UNIT=io_unit, FMT="(/,T2,A,/,/,T3,A,I0,/,/,T3,A2,I0,A4,F10.6)") &
     883        2090 :                   "Parameters for the always stable predictor-corrector (ASPC) method:", &
     884        2090 :                   "ASPC order: ", MAX(nvec - 2, 0), &
     885        4180 :                   "B(", 1, ") = ", alpha
     886             :             END IF
     887       10921 :             CALL cp_fm_scale(alpha, mo_coeff)
     888             : 
     889       41411 :             DO i = 2, nvec
     890       30490 :                t0_state => wfi_get_snapshot(wf_history, wf_index=i)
     891       30490 :                IF (use_overlap) THEN
     892       18722 :                   CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
     893       18722 :                   CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
     894             :                ELSE
     895             :                   CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
     896       11768 :                                      t1_state%wf(ispin), 0.0_dp, csc)
     897             :                END IF
     898       30490 :                CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
     899             :                alpha = (-1.0_dp)**(i + 1)*REAL(i, KIND=dp)* &
     900       30490 :                        binomial(2*nvec, nvec - i)/binomial(2*nvec - 2, nvec - 1)
     901       30490 :                IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
     902             :                   WRITE (UNIT=io_unit, FMT="(T3,A2,I0,A4,F10.6)") &
     903        5878 :                      "B(", i, ") = ", alpha
     904             :                END IF
     905       41411 :                CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp)
     906             :             END DO
     907       10921 :             CALL cp_fm_release(csc)
     908       10921 :             CALL cp_fm_release(fm_tmp)
     909             :             CALL reorthogonalize_vectors(qs_env, &
     910             :                                          v_matrix=mo_coeff, &
     911       10921 :                                          n_col=k)
     912             :             CALL calculate_density_matrix(mo_set=mos(ispin), &
     913       32302 :                                           density_matrix=rho_ao(ispin)%matrix)
     914             :          END DO
     915       10460 :          CALL qs_rho_update_rho(rho, qs_env=qs_env)
     916       10460 :          CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
     917             : 
     918             :       CASE default
     919             :          CALL cp_abort(__LOCATION__, &
     920             :                        "Unknown interpolation method: "// &
     921       17495 :                        TRIM(ADJUSTL(cp_to_string(wf_history%interpolation_method_nr))))
     922             :       END SELECT
     923       17495 :       IF (PRESENT(orthogonal_wf)) orthogonal_wf = my_orthogonal_wf
     924             :       CALL cp_print_key_finished_output(io_unit, logger, qs_env%input, &
     925       17495 :                                         "DFT%SCF%PRINT%PROGRAM_RUN_INFO")
     926       17495 :       CALL timestop(handle)
     927       17495 :    END SUBROUTINE wfi_extrapolate
     928             : 
     929             : ! **************************************************************************************************
     930             : !> \brief Decides if scf control variables has to changed due
     931             : !>      to using a WF extrapolation.
     932             : !> \param qs_env The QS environment
     933             : !> \param nvec ...
     934             : !> \par History
     935             : !>      11.2006 created [TdK]
     936             : !> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch)
     937             : ! **************************************************************************************************
     938        6775 :    ELEMENTAL SUBROUTINE wfi_set_history_variables(qs_env, nvec)
     939             :       TYPE(qs_environment_type), INTENT(INOUT)           :: qs_env
     940             :       INTEGER, INTENT(IN)                                :: nvec
     941             : 
     942        6775 :       IF (nvec >= qs_env%wf_history%memory_depth) THEN
     943        1547 :          IF ((qs_env%scf_control%max_scf_hist .NE. 0) .AND. (qs_env%scf_control%eps_scf_hist .NE. 0)) THEN
     944           0 :             qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
     945           0 :             qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
     946           0 :             qs_env%scf_control%outer_scf%have_scf = .FALSE.
     947        1547 :          ELSE IF (qs_env%scf_control%max_scf_hist .NE. 0) THEN
     948           0 :             qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
     949           0 :             qs_env%scf_control%outer_scf%have_scf = .FALSE.
     950        1547 :          ELSE IF (qs_env%scf_control%eps_scf_hist .NE. 0) THEN
     951           0 :             qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
     952           0 :             qs_env%scf_control%outer_scf%eps_scf = qs_env%scf_control%eps_scf_hist
     953             :          END IF
     954             :       END IF
     955             : 
     956        6775 :    END SUBROUTINE wfi_set_history_variables
     957             : 
     958             : ! **************************************************************************************************
     959             : !> \brief updates the snapshot buffer, taking a new snapshot
     960             : !> \param wf_history the history buffer to update
     961             : !> \param qs_env the qs_env we get the info from
     962             : !> \param dt ...
     963             : !> \par History
     964             : !>      02.2003 created [fawzi]
     965             : !> \author fawzi
     966             : ! **************************************************************************************************
     967       17499 :    SUBROUTINE wfi_update(wf_history, qs_env, dt)
     968             :       TYPE(qs_wf_history_type), POINTER                  :: wf_history
     969             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     970             :       REAL(KIND=dp), INTENT(in)                          :: dt
     971             : 
     972       17499 :       CPASSERT(ASSOCIATED(wf_history))
     973       17499 :       CPASSERT(wf_history%ref_count > 0)
     974       17499 :       CPASSERT(ASSOCIATED(qs_env))
     975             : 
     976       17499 :       wf_history%snapshot_count = wf_history%snapshot_count + 1
     977       17499 :       IF (wf_history%memory_depth > 0) THEN
     978             :          wf_history%last_state_index = MODULO(wf_history%snapshot_count, &
     979       16450 :                                               wf_history%memory_depth) + 1
     980             :          CALL wfs_update(snapshot=wf_history%past_states &
     981             :                          (wf_history%last_state_index)%snapshot, wf_history=wf_history, &
     982       16450 :                          qs_env=qs_env, dt=dt)
     983             :       END IF
     984       17499 :    END SUBROUTINE wfi_update
     985             : 
     986             : ! **************************************************************************************************
     987             : !> \brief reorthogonalizes the mos
     988             : !> \param qs_env the qs_env in which to orthogonalize
     989             : !> \param v_matrix the vectors to orthogonalize
     990             : !> \param n_col number of column of v to orthogonalize
     991             : !> \par History
     992             : !>      04.2003 created [fawzi]
     993             : !> \author Fawzi Mohamed
     994             : ! **************************************************************************************************
     995       23462 :    SUBROUTINE reorthogonalize_vectors(qs_env, v_matrix, n_col)
     996             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     997             :       TYPE(cp_fm_type), INTENT(IN)                       :: v_matrix
     998             :       INTEGER, INTENT(in), OPTIONAL                      :: n_col
     999             : 
    1000             :       CHARACTER(len=*), PARAMETER :: routineN = 'reorthogonalize_vectors'
    1001             : 
    1002             :       INTEGER                                            :: handle, my_n_col
    1003             :       LOGICAL                                            :: has_unit_metric, &
    1004             :                                                             ortho_contains_cholesky, &
    1005             :                                                             smearing_is_used
    1006             :       TYPE(cp_fm_pool_type), POINTER                     :: maxao_maxmo_fm_pool
    1007       11731 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
    1008             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1009             :       TYPE(qs_matrix_pools_type), POINTER                :: mpools
    1010             :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
    1011             :       TYPE(scf_control_type), POINTER                    :: scf_control
    1012             : 
    1013       11731 :       NULLIFY (scf_env, scf_control, maxao_maxmo_fm_pool, matrix_s, mpools, dft_control)
    1014       11731 :       CALL timeset(routineN, handle)
    1015             : 
    1016       11731 :       CPASSERT(ASSOCIATED(qs_env))
    1017             : 
    1018       11731 :       CALL cp_fm_get_info(v_matrix, ncol_global=my_n_col)
    1019       11731 :       IF (PRESENT(n_col)) my_n_col = n_col
    1020             :       CALL get_qs_env(qs_env, mpools=mpools, &
    1021             :                       scf_env=scf_env, &
    1022             :                       scf_control=scf_control, &
    1023             :                       matrix_s=matrix_s, &
    1024       11731 :                       dft_control=dft_control)
    1025       11731 :       CALL mpools_get(mpools, maxao_maxmo_fm_pool=maxao_maxmo_fm_pool)
    1026       11731 :       IF (ASSOCIATED(scf_env)) THEN
    1027             :          ortho_contains_cholesky = (scf_env%method /= ot_method_nr) .AND. &
    1028             :                                    (scf_env%cholesky_method > 0) .AND. &
    1029       11731 :                                    ASSOCIATED(scf_env%ortho)
    1030             :       ELSE
    1031             :          ortho_contains_cholesky = .FALSE.
    1032             :       END IF
    1033             : 
    1034       11731 :       CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
    1035       11731 :       smearing_is_used = .FALSE.
    1036       11731 :       IF (dft_control%smear) THEN
    1037          94 :          smearing_is_used = .TRUE.
    1038             :       END IF
    1039             : 
    1040       11731 :       IF (has_unit_metric) THEN
    1041        3462 :          CALL make_basis_simple(v_matrix, my_n_col)
    1042        8269 :       ELSE IF (smearing_is_used) THEN
    1043             :          CALL make_basis_lowdin(vmatrix=v_matrix, ncol=my_n_col, &
    1044          94 :                                 matrix_s=matrix_s(1)%matrix)
    1045        8175 :       ELSE IF (ortho_contains_cholesky) THEN
    1046             :          CALL make_basis_cholesky(vmatrix=v_matrix, ncol=my_n_col, &
    1047        5264 :                                   ortho=scf_env%ortho)
    1048             :       ELSE
    1049        2911 :          CALL make_basis_sm(v_matrix, my_n_col, matrix_s(1)%matrix)
    1050             :       END IF
    1051       11731 :       CALL timestop(handle)
    1052       11731 :    END SUBROUTINE reorthogonalize_vectors
    1053             : 
    1054             : ! **************************************************************************************************
    1055             : !> \brief purges wf_history retaining only the latest snapshot
    1056             : !> \param qs_env the qs env with the latest result, and that will contain
    1057             : !>        the purged wf_history
    1058             : !> \par History
    1059             : !>      05.2016 created [Nico Holmberg]
    1060             : !> \author Nico Holmberg
    1061             : ! **************************************************************************************************
    1062           0 :    SUBROUTINE wfi_purge_history(qs_env)
    1063             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1064             : 
    1065             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'wfi_purge_history'
    1066             : 
    1067             :       INTEGER                                            :: handle, io_unit, print_level
    1068             :       TYPE(cp_logger_type), POINTER                      :: logger
    1069             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1070             :       TYPE(qs_wf_history_type), POINTER                  :: wf_history
    1071             : 
    1072           0 :       NULLIFY (dft_control, wf_history)
    1073             : 
    1074           0 :       CALL timeset(routineN, handle)
    1075           0 :       logger => cp_get_default_logger()
    1076           0 :       print_level = logger%iter_info%print_level
    1077             :       io_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%SCF%PRINT%PROGRAM_RUN_INFO", &
    1078           0 :                                      extension=".scfLog")
    1079             : 
    1080           0 :       CPASSERT(ASSOCIATED(qs_env))
    1081           0 :       CPASSERT(ASSOCIATED(qs_env%wf_history))
    1082           0 :       CPASSERT(qs_env%wf_history%ref_count > 0)
    1083           0 :       CALL get_qs_env(qs_env, dft_control=dft_control)
    1084             : 
    1085           0 :       SELECT CASE (qs_env%wf_history%interpolation_method_nr)
    1086             :       CASE (wfi_use_guess_method_nr, wfi_use_prev_wf_method_nr, &
    1087             :             wfi_use_prev_p_method_nr, wfi_use_prev_rho_r_method_nr, &
    1088             :             wfi_frozen_method_nr)
    1089             :          ! do nothing
    1090             :       CASE (wfi_linear_wf_method_nr, wfi_linear_p_method_nr, &
    1091             :             wfi_linear_ps_method_nr, wfi_ps_method_nr, &
    1092             :             wfi_aspc_nr)
    1093           0 :          IF (qs_env%wf_history%snapshot_count .GE. 2) THEN
    1094           0 :             IF (debug_this_module .AND. io_unit > 0) &
    1095           0 :                WRITE (io_unit, FMT="(T2,A)") "QS| Purging WFN history"
    1096             :             CALL wfi_create(wf_history, interpolation_method_nr= &
    1097             :                             dft_control%qs_control%wf_interpolation_method_nr, &
    1098             :                             extrapolation_order=dft_control%qs_control%wf_extrapolation_order, &
    1099           0 :                             has_unit_metric=qs_env%has_unit_metric)
    1100             :             CALL set_qs_env(qs_env=qs_env, &
    1101           0 :                             wf_history=wf_history)
    1102           0 :             CALL wfi_release(wf_history)
    1103           0 :             CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
    1104             :          END IF
    1105             :       CASE DEFAULT
    1106           0 :          CPABORT("Unknown extrapolation method.")
    1107             :       END SELECT
    1108           0 :       CALL timestop(handle)
    1109             : 
    1110           0 :    END SUBROUTINE wfi_purge_history
    1111             : 
    1112             : END MODULE qs_wf_history_methods

Generated by: LCOV version 1.15