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 cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale,&
27 : cp_cfm_gemm,&
28 : cp_cfm_scale_and_add,&
29 : cp_cfm_scale_and_add_fm,&
30 : cp_cfm_triangular_multiply
31 : USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose
32 : USE cp_cfm_diag, ONLY: cp_cfm_heevd
33 : USE cp_cfm_types, ONLY: cp_cfm_create,&
34 : cp_cfm_release,&
35 : cp_cfm_set_all,&
36 : cp_cfm_to_cfm,&
37 : cp_cfm_to_fm,&
38 : cp_cfm_type,&
39 : 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_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_release, cp_fm_set_all, cp_fm_start_copy_general, cp_fm_to_fm, &
66 : 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_rho_r_method_nr, wfi_use_prev_wf_method_nr
78 : USE kinds, ONLY: dp
79 : USE kpoint_methods, ONLY: kpoint_density_matrices,&
80 : kpoint_density_transform,&
81 : kpoint_set_mo_occupation,&
82 : rskp_transform
83 : USE kpoint_types, ONLY: get_kpoint_info,&
84 : kpoint_env_type,&
85 : kpoint_type
86 : USE mathconstants, ONLY: gaussi,&
87 : z_one,&
88 : z_zero
89 : USE mathlib, ONLY: binomial
90 : USE message_passing, ONLY: mp_para_env_type
91 : USE parallel_gemm_api, ONLY: parallel_gemm
92 : USE pw_env_types, ONLY: pw_env_get,&
93 : pw_env_type
94 : USE pw_methods, ONLY: pw_copy
95 : USE pw_pool_types, ONLY: pw_pool_type
96 : USE pw_types, ONLY: pw_c1d_gs_type,&
97 : pw_r3d_rs_type
98 : USE qs_density_matrices, ONLY: calculate_density_matrix
99 : USE qs_environment_types, ONLY: get_qs_env,&
100 : qs_environment_type,&
101 : set_qs_env
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_set,&
115 : qs_rho_type
116 : USE qs_scf_types, ONLY: ot_method_nr,&
117 : qs_scf_env_type
118 : USE qs_wf_history_types, ONLY: qs_wf_history_type,&
119 : qs_wf_snapshot_type,&
120 : wfi_get_snapshot,&
121 : wfi_release
122 : USE scf_control_types, ONLY: scf_control_type
123 : #include "./base/base_uses.f90"
124 :
125 : IMPLICIT NONE
126 : PRIVATE
127 :
128 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
129 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_wf_history_methods'
130 :
131 : PUBLIC :: wfi_create, wfi_update, wfi_create_for_kp, &
132 : wfi_extrapolate, wfi_get_method_label, &
133 : reorthogonalize_vectors, wfi_purge_history
134 :
135 : CONTAINS
136 :
137 : ! **************************************************************************************************
138 : !> \brief allocates and initialize a wavefunction snapshot
139 : !> \param snapshot the snapshot to create
140 : !> \par History
141 : !> 02.2003 created [fawzi]
142 : !> 02.2005 added wf_mol [MI]
143 : !> \author fawzi
144 : ! **************************************************************************************************
145 11434 : SUBROUTINE wfs_create(snapshot)
146 : TYPE(qs_wf_snapshot_type), INTENT(OUT) :: snapshot
147 :
148 : NULLIFY (snapshot%wf, snapshot%rho_r, &
149 : snapshot%rho_g, snapshot%rho_ao, snapshot%rho_ao_kp, &
150 : snapshot%overlap, snapshot%wf_kp, snapshot%overlap_cfm_kp, &
151 : snapshot%rho_frozen)
152 11434 : snapshot%dt = 1.0_dp
153 11434 : END SUBROUTINE wfs_create
154 :
155 : ! **************************************************************************************************
156 : !> \brief updates the given snapshot
157 : !> \param snapshot the snapshot to be updated
158 : !> \param wf_history the history
159 : !> \param qs_env the qs_env that should be snapshotted
160 : !> \param dt the time of the snapshot (wrt. to the previous snapshot)
161 : !> \par History
162 : !> 02.2003 created [fawzi]
163 : !> 02.2005 added kg_fm_mol_set for KG_GPW [MI]
164 : !> \author fawzi
165 : ! **************************************************************************************************
166 20512 : SUBROUTINE wfs_update(snapshot, wf_history, qs_env, dt)
167 : TYPE(qs_wf_snapshot_type), POINTER :: snapshot
168 : TYPE(qs_wf_history_type), POINTER :: wf_history
169 : TYPE(qs_environment_type), POINTER :: qs_env
170 : REAL(KIND=dp), INTENT(in), OPTIONAL :: dt
171 :
172 : CHARACTER(len=*), PARAMETER :: routineN = 'wfs_update'
173 :
174 : INTEGER :: handle, ic, igroup, ik, ikp, img, &
175 : indx_ft, ispin, kplocal, nc, nimg, &
176 : nkp_all, nkp_grps, nspin_kp, nspins
177 : INTEGER, DIMENSION(2) :: kp_range
178 20512 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
179 20512 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
180 : LOGICAL :: my_kpgrp
181 20512 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
182 20512 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info_ft
183 20512 : 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 20512 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
188 20512 : 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 20512 : 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 20512 : POINTER :: sab_nl
197 20512 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
198 : TYPE(pw_env_type), POINTER :: pw_env
199 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
200 20512 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
201 : TYPE(qs_matrix_pools_type), POINTER :: mpools_kp
202 : TYPE(qs_rho_type), POINTER :: rho
203 : TYPE(qs_scf_env_type), POINTER :: scf_env
204 :
205 20512 : CALL timeset(routineN, handle)
206 :
207 20512 : NULLIFY (pw_env, auxbas_pw_pool, ao_mo_pools, ao_ao_fm_pools, dft_control, mos, mo_coeff, &
208 20512 : rho, rho_r, rho_g, rho_ao, matrix_s, matrix_s_kp, kpoints, kp, &
209 20512 : kp_dist, cell_to_index, xkp, sab_nl, scf_env, mpools_kp, para_env_ft, &
210 20512 : rmat_ft, cmat_ft, tmpmat_ft, ao_ao_struct_ft)
211 : CALL get_qs_env(qs_env, pw_env=pw_env, &
212 20512 : dft_control=dft_control, rho=rho)
213 20512 : CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_pools)
214 20512 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
215 :
216 20512 : CPASSERT(ASSOCIATED(wf_history))
217 20512 : CPASSERT(ASSOCIATED(dft_control))
218 20512 : IF (.NOT. ASSOCIATED(snapshot)) THEN
219 11434 : ALLOCATE (snapshot)
220 11434 : CALL wfs_create(snapshot)
221 : END IF
222 20512 : CPASSERT(wf_history%ref_count > 0)
223 :
224 20512 : nspins = dft_control%nspins
225 20512 : snapshot%dt = 1.0_dp
226 20512 : IF (PRESENT(dt)) snapshot%dt = dt
227 20512 : IF (wf_history%store_wf) THEN
228 19462 : CALL get_qs_env(qs_env, mos=mos)
229 19462 : IF (.NOT. ASSOCIATED(snapshot%wf)) THEN
230 : CALL fm_pools_create_fm_vect(ao_mo_pools, snapshot%wf, &
231 11152 : name="ws_snap-ws")
232 11152 : CPASSERT(nspins == SIZE(snapshot%wf))
233 : END IF
234 41066 : DO ispin = 1, nspins
235 21604 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
236 41066 : CALL cp_fm_to_fm(mo_coeff, snapshot%wf(ispin))
237 : END DO
238 : ELSE
239 1050 : CALL fm_pools_give_back_fm_vect(ao_mo_pools, snapshot%wf)
240 : END IF
241 :
242 20512 : IF (wf_history%store_rho_r) THEN
243 0 : CALL qs_rho_get(rho, rho_r=rho_r)
244 0 : CPASSERT(ASSOCIATED(rho_r))
245 0 : IF (.NOT. ASSOCIATED(snapshot%rho_r)) THEN
246 0 : ALLOCATE (snapshot%rho_r(nspins))
247 0 : DO ispin = 1, nspins
248 0 : CALL auxbas_pw_pool%create_pw(snapshot%rho_r(ispin))
249 : END DO
250 : END IF
251 0 : DO ispin = 1, nspins
252 0 : CALL pw_copy(rho_r(ispin), snapshot%rho_r(ispin))
253 : END DO
254 20512 : ELSE IF (ASSOCIATED(snapshot%rho_r)) THEN
255 0 : DO ispin = 1, SIZE(snapshot%rho_r)
256 0 : CALL auxbas_pw_pool%give_back_pw(snapshot%rho_r(ispin))
257 : END DO
258 0 : DEALLOCATE (snapshot%rho_r)
259 : END IF
260 :
261 20512 : IF (wf_history%store_rho_g) THEN
262 0 : CALL qs_rho_get(rho, rho_g=rho_g)
263 0 : CPASSERT(ASSOCIATED(rho_g))
264 0 : IF (.NOT. ASSOCIATED(snapshot%rho_g)) THEN
265 0 : ALLOCATE (snapshot%rho_g(nspins))
266 0 : DO ispin = 1, nspins
267 0 : CALL auxbas_pw_pool%create_pw(snapshot%rho_g(ispin))
268 : END DO
269 : END IF
270 0 : DO ispin = 1, nspins
271 0 : CALL pw_copy(rho_g(ispin), snapshot%rho_g(ispin))
272 : END DO
273 20512 : ELSE IF (ASSOCIATED(snapshot%rho_g)) THEN
274 0 : DO ispin = 1, SIZE(snapshot%rho_g)
275 0 : CALL auxbas_pw_pool%give_back_pw(snapshot%rho_g(ispin))
276 : END DO
277 0 : DEALLOCATE (snapshot%rho_g)
278 : END IF
279 :
280 20512 : IF (ASSOCIATED(snapshot%rho_ao)) THEN ! the sparsity might be different
281 : ! (future struct:check)
282 270 : CALL dbcsr_deallocate_matrix_set(snapshot%rho_ao)
283 : END IF
284 20512 : IF (wf_history%store_rho_ao) THEN
285 338 : CALL qs_rho_get(rho, rho_ao=rho_ao)
286 338 : CPASSERT(ASSOCIATED(rho_ao))
287 :
288 338 : CALL dbcsr_allocate_matrix_set(snapshot%rho_ao, nspins)
289 836 : DO ispin = 1, nspins
290 498 : ALLOCATE (snapshot%rho_ao(ispin)%matrix)
291 836 : CALL dbcsr_copy(snapshot%rho_ao(ispin)%matrix, rho_ao(ispin)%matrix)
292 : END DO
293 : END IF
294 :
295 20512 : IF (ASSOCIATED(snapshot%rho_ao_kp)) THEN ! the sparsity might be different
296 : ! (future struct:check)
297 220 : CALL dbcsr_deallocate_matrix_set(snapshot%rho_ao_kp)
298 : END IF
299 20512 : IF (wf_history%store_rho_ao_kp) THEN
300 232 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
301 232 : CPASSERT(ASSOCIATED(rho_ao_kp))
302 :
303 232 : nimg = dft_control%nimages
304 232 : CALL dbcsr_allocate_matrix_set(snapshot%rho_ao_kp, nspins, nimg)
305 554 : DO ispin = 1, nspins
306 34092 : DO img = 1, nimg
307 33538 : ALLOCATE (snapshot%rho_ao_kp(ispin, img)%matrix)
308 : CALL dbcsr_copy(snapshot%rho_ao_kp(ispin, img)%matrix, &
309 33860 : rho_ao_kp(ispin, img)%matrix)
310 : END DO
311 : END DO
312 : END IF
313 :
314 20512 : IF (ASSOCIATED(snapshot%overlap)) THEN ! the sparsity might be different
315 : ! (future struct:check)
316 5702 : CALL dbcsr_deallocate_matrix(snapshot%overlap)
317 : END IF
318 20512 : IF (wf_history%store_overlap) THEN
319 15820 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
320 15820 : CPASSERT(ASSOCIATED(matrix_s))
321 15820 : CPASSERT(ASSOCIATED(matrix_s(1)%matrix))
322 15820 : ALLOCATE (snapshot%overlap)
323 15820 : CALL dbcsr_copy(snapshot%overlap, matrix_s(1)%matrix)
324 : END IF
325 :
326 20512 : CALL get_qs_env(qs_env, kpoints=kpoints)
327 20512 : IF (ASSOCIATED(kpoints)) THEN
328 20512 : IF (ASSOCIATED(kpoints%kp_env)) THEN
329 : ! --- k-point WFN snapshot: store complex MO coefficients per local k-point ---
330 696 : IF (wf_history%store_wf_kp) THEN
331 464 : CALL get_kpoint_info(kpoints, kp_range=kp_range)
332 464 : kplocal = kp_range(2) - kp_range(1) + 1
333 464 : nspin_kp = SIZE(kpoints%kp_env(1)%kpoint_env%mos, 2)
334 464 : nc = SIZE(kpoints%kp_env(1)%kpoint_env%mos, 1) ! 2=complex, 1=real
335 :
336 464 : IF (ASSOCIATED(snapshot%wf_kp)) THEN
337 732 : DO ikp = 1, SIZE(snapshot%wf_kp, 1)
338 1664 : DO ic = 1, SIZE(snapshot%wf_kp, 2)
339 2330 : DO ispin = 1, SIZE(snapshot%wf_kp, 3)
340 1864 : CALL cp_fm_release(snapshot%wf_kp(ikp, ic, ispin))
341 : END DO
342 : END DO
343 : END DO
344 266 : DEALLOCATE (snapshot%wf_kp)
345 : END IF
346 :
347 7474 : ALLOCATE (snapshot%wf_kp(kplocal, nc, nspin_kp))
348 1970 : DO ikp = 1, kplocal
349 1506 : kp => kpoints%kp_env(ikp)%kpoint_env
350 3824 : DO ispin = 1, nspin_kp
351 7068 : DO ic = 1, nc
352 3708 : CALL get_mo_set(kp%mos(ic, ispin), mo_coeff=mo_coeff)
353 : CALL cp_fm_create(snapshot%wf_kp(ikp, ic, ispin), &
354 : mo_coeff%matrix_struct, &
355 3708 : name="wfkp_snap")
356 5562 : CALL cp_fm_to_fm(mo_coeff, snapshot%wf_kp(ikp, ic, ispin))
357 : END DO
358 : END DO
359 : END DO
360 : END IF
361 :
362 : ! --- k-point overlap snapshot: Fourier-transform S(R)→S(k) NOW and store as cfm ---
363 : ! This is critical: we MUST transform at snapshot time using the CURRENT neighbor
364 : ! list. Storing S(R) and re-transforming later would use a stale neighbor list,
365 : ! producing wrong S(k) when the neighbor list changes during MD.
366 696 : IF (wf_history%store_overlap_kp) THEN
367 464 : CALL get_qs_env(qs_env, matrix_s_kp=matrix_s_kp, scf_env=scf_env)
368 : CALL get_kpoint_info(kpoints, nkp=nkp_all, xkp=xkp, kp_range=kp_range, &
369 : nkp_groups=nkp_grps, kp_dist=kp_dist, &
370 464 : sab_nl=sab_nl, cell_to_index=cell_to_index)
371 464 : kplocal = kp_range(2) - kp_range(1) + 1
372 464 : para_env_ft => kpoints%blacs_env_all%para_env
373 :
374 : ! Allocate dbcsr work matrices for FT (same pattern as do_general_diag_kp)
375 464 : ALLOCATE (rmat_ft, cmat_ft, tmpmat_ft)
376 : CALL dbcsr_create(rmat_ft, template=matrix_s_kp(1, 1)%matrix, &
377 464 : matrix_type=dbcsr_type_symmetric)
378 : CALL dbcsr_create(cmat_ft, template=matrix_s_kp(1, 1)%matrix, &
379 464 : matrix_type=dbcsr_type_antisymmetric)
380 : CALL dbcsr_create(tmpmat_ft, template=matrix_s_kp(1, 1)%matrix, &
381 464 : matrix_type=dbcsr_type_no_symmetry)
382 464 : CALL cp_dbcsr_alloc_block_from_nbl(rmat_ft, sab_nl)
383 464 : CALL cp_dbcsr_alloc_block_from_nbl(cmat_ft, sab_nl)
384 :
385 : ! Get kp-subgroup FM from pool
386 464 : CALL get_kpoint_info(kpoints, mpools=mpools_kp)
387 464 : CALL mpools_get(mpools_kp, ao_ao_fm_pools=ao_ao_fm_pools)
388 464 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal_ft)
389 :
390 : ! Release old snapshot if present
391 464 : IF (ASSOCIATED(snapshot%overlap_cfm_kp)) THEN
392 732 : DO ikp = 1, SIZE(snapshot%overlap_cfm_kp)
393 732 : CALL cp_cfm_release(snapshot%overlap_cfm_kp(ikp))
394 : END DO
395 266 : DEALLOCATE (snapshot%overlap_cfm_kp)
396 : END IF
397 2898 : ALLOCATE (snapshot%overlap_cfm_kp(kplocal))
398 :
399 464 : CALL cp_fm_get_info(fmlocal_ft, matrix_struct=ao_ao_struct_ft)
400 :
401 : ! Communication info array
402 11112 : ALLOCATE (info_ft(kplocal*nkp_grps, 2))
403 :
404 : ! Phase A: Start async FT + redistribute for each k-point
405 464 : indx_ft = 0
406 1970 : DO ikp = 1, kplocal
407 4510 : DO igroup = 1, nkp_grps
408 2540 : ik = kp_dist(1, igroup) + ikp - 1
409 2540 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
410 2540 : indx_ft = indx_ft + 1
411 :
412 2540 : CALL dbcsr_set(rmat_ft, 0.0_dp)
413 2540 : CALL dbcsr_set(cmat_ft, 0.0_dp)
414 : CALL rskp_transform(rmatrix=rmat_ft, cmatrix=cmat_ft, rsmat=matrix_s_kp, &
415 : ispin=1, xkp=xkp(1:3, ik), &
416 2540 : cell_to_index=cell_to_index, sab_nl=sab_nl)
417 2540 : CALL dbcsr_desymmetrize(rmat_ft, tmpmat_ft)
418 2540 : CALL copy_dbcsr_to_fm(tmpmat_ft, scf_env%scf_work1(1))
419 2540 : CALL dbcsr_desymmetrize(cmat_ft, tmpmat_ft)
420 2540 : CALL copy_dbcsr_to_fm(tmpmat_ft, scf_env%scf_work1(2))
421 :
422 4046 : IF (my_kpgrp) THEN
423 : CALL cp_fm_start_copy_general(scf_env%scf_work1(1), fmlocal_ft, &
424 1506 : para_env_ft, info_ft(indx_ft, 1))
425 : CALL cp_fm_start_copy_general(scf_env%scf_work1(2), fmlocal_ft, &
426 1506 : para_env_ft, info_ft(indx_ft, 2))
427 : ELSE
428 : CALL cp_fm_start_copy_general(scf_env%scf_work1(1), fmdummy_ft, &
429 1034 : para_env_ft, info_ft(indx_ft, 1))
430 : CALL cp_fm_start_copy_general(scf_env%scf_work1(2), fmdummy_ft, &
431 1034 : para_env_ft, info_ft(indx_ft, 2))
432 : END IF
433 : END DO
434 : END DO
435 :
436 : ! Phase B: Finish communication and assemble S(k) as cfm
437 : indx_ft = 0
438 1970 : DO ikp = 1, kplocal
439 1506 : CALL cp_cfm_create(snapshot%overlap_cfm_kp(ikp), ao_ao_struct_ft)
440 1506 : CALL cp_cfm_set_all(snapshot%overlap_cfm_kp(ikp), z_zero)
441 4510 : DO igroup = 1, nkp_grps
442 2540 : ik = kp_dist(1, igroup) + ikp - 1
443 2540 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
444 1034 : indx_ft = indx_ft + 1
445 1506 : IF (my_kpgrp) THEN
446 1506 : CALL cp_fm_finish_copy_general(fmlocal_ft, info_ft(indx_ft, 1))
447 : CALL cp_cfm_scale_and_add_fm(z_zero, snapshot%overlap_cfm_kp(ikp), &
448 1506 : z_one, fmlocal_ft)
449 1506 : CALL cp_fm_finish_copy_general(fmlocal_ft, info_ft(indx_ft, 2))
450 : CALL cp_cfm_scale_and_add_fm(z_one, snapshot%overlap_cfm_kp(ikp), &
451 1506 : gaussi, fmlocal_ft)
452 : END IF
453 : END DO
454 : END DO
455 :
456 : ! Cleanup
457 3004 : DO indx_ft = 1, kplocal*nkp_grps
458 2540 : CALL cp_fm_cleanup_copy_general(info_ft(indx_ft, 1))
459 3004 : CALL cp_fm_cleanup_copy_general(info_ft(indx_ft, 2))
460 : END DO
461 6008 : DEALLOCATE (info_ft)
462 464 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal_ft)
463 464 : CALL dbcsr_deallocate_matrix(rmat_ft)
464 464 : CALL dbcsr_deallocate_matrix(cmat_ft)
465 928 : CALL dbcsr_deallocate_matrix(tmpmat_ft)
466 : END IF
467 : END IF
468 : END IF
469 :
470 : IF (wf_history%store_frozen_density) THEN
471 : ! do nothing
472 : ! CALL deallocate_matrix_set(snapshot%rho_frozen%rho_ao)
473 : END IF
474 :
475 20512 : CALL timestop(handle)
476 :
477 20512 : END SUBROUTINE wfs_update
478 :
479 : ! **************************************************************************************************
480 : !> \brief ...
481 : !> \param wf_history ...
482 : !> \param interpolation_method_nr the tag of the method used for
483 : !> the extrapolation of the initial density for the next md step
484 : !> (see qs_wf_history_types:wfi_*_method_nr)
485 : !> \param extrapolation_order ...
486 : !> \param has_unit_metric ...
487 : !> \par History
488 : !> 02.2003 created [fawzi]
489 : !> \author fawzi
490 : ! **************************************************************************************************
491 7748 : SUBROUTINE wfi_create(wf_history, interpolation_method_nr, extrapolation_order, &
492 : has_unit_metric)
493 : TYPE(qs_wf_history_type), POINTER :: wf_history
494 : INTEGER, INTENT(in) :: interpolation_method_nr, &
495 : extrapolation_order
496 : LOGICAL, INTENT(IN) :: has_unit_metric
497 :
498 : CHARACTER(len=*), PARAMETER :: routineN = 'wfi_create'
499 :
500 : INTEGER :: handle, i
501 :
502 7748 : CALL timeset(routineN, handle)
503 :
504 7748 : ALLOCATE (wf_history)
505 7748 : wf_history%ref_count = 1
506 7748 : wf_history%memory_depth = 0
507 7748 : wf_history%snapshot_count = 0
508 7748 : wf_history%last_state_index = 1
509 : wf_history%store_wf = .FALSE.
510 : wf_history%store_rho_r = .FALSE.
511 : wf_history%store_rho_g = .FALSE.
512 : wf_history%store_rho_ao = .FALSE.
513 : wf_history%store_rho_ao_kp = .FALSE.
514 : wf_history%store_overlap = .FALSE.
515 : wf_history%store_wf_kp = .FALSE.
516 : wf_history%store_overlap_kp = .FALSE.
517 : wf_history%store_frozen_density = .FALSE.
518 : NULLIFY (wf_history%past_states)
519 :
520 7748 : wf_history%interpolation_method_nr = interpolation_method_nr
521 :
522 : SELECT CASE (wf_history%interpolation_method_nr)
523 : CASE (wfi_use_guess_method_nr)
524 : wf_history%memory_depth = 0
525 : CASE (wfi_use_prev_wf_method_nr)
526 64 : wf_history%memory_depth = 0
527 : CASE (wfi_use_prev_p_method_nr)
528 64 : wf_history%memory_depth = 1
529 64 : wf_history%store_rho_ao = .TRUE.
530 : CASE (wfi_use_prev_rho_r_method_nr)
531 4 : wf_history%memory_depth = 1
532 4 : 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 343 : CALL cite_reference(VandeVondele2005a)
545 343 : wf_history%memory_depth = extrapolation_order + 1
546 343 : wf_history%store_wf = .TRUE.
547 343 : wf_history%store_wf_kp = .TRUE.
548 343 : IF (.NOT. has_unit_metric) THEN
549 339 : wf_history%store_overlap = .TRUE.
550 339 : 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 7171 : wf_history%memory_depth = extrapolation_order + 2
557 7171 : wf_history%store_wf = .TRUE.
558 7171 : wf_history%store_wf_kp = .TRUE.
559 7171 : IF (.NOT. has_unit_metric) THEN
560 6189 : wf_history%store_overlap = .TRUE.
561 6189 : wf_history%store_overlap_kp = .TRUE.
562 : END IF
563 : CASE (wfi_gext_proj_nr)
564 4 : wf_history%memory_depth = extrapolation_order
565 4 : wf_history%store_wf = .TRUE.
566 4 : wf_history%store_overlap = .TRUE.
567 : CASE (wfi_gext_proj_qtr_nr)
568 4 : wf_history%memory_depth = extrapolation_order
569 4 : wf_history%store_wf = .TRUE.
570 4 : wf_history%store_overlap = .TRUE.
571 : CASE default
572 : CALL cp_abort(__LOCATION__, &
573 : "Unknown interpolation method: "// &
574 0 : TRIM(ADJUSTL(cp_to_string(interpolation_method_nr))))
575 7748 : wf_history%interpolation_method_nr = wfi_use_prev_rho_r_method_nr
576 : END SELECT
577 60157 : ALLOCATE (wf_history%past_states(wf_history%memory_depth))
578 :
579 44799 : DO i = 1, SIZE(wf_history%past_states)
580 44799 : NULLIFY (wf_history%past_states(i)%snapshot)
581 : END DO
582 :
583 7748 : CALL timestop(handle)
584 7748 : END SUBROUTINE wfi_create
585 :
586 : ! **************************************************************************************************
587 : !> \brief Adapts wf_history storage flags for k-point calculations.
588 : !> For ASPC, switches from Gamma WFN storage to k-point WFN storage.
589 : !> Other WFN-based methods remain blocked.
590 : !> \param wf_history ...
591 : !> \par History
592 : !> 06.2015 created [jhu]
593 : !> \author jhu
594 : ! **************************************************************************************************
595 198 : SUBROUTINE wfi_create_for_kp(wf_history)
596 : TYPE(qs_wf_history_type), POINTER :: wf_history
597 :
598 198 : CPASSERT(ASSOCIATED(wf_history))
599 198 : IF (wf_history%store_rho_ao) THEN
600 10 : wf_history%store_rho_ao_kp = .TRUE.
601 10 : wf_history%store_rho_ao = .FALSE.
602 : END IF
603 : ! KP-compatible PS and ASPC: switch from Gamma WFN storage to KP WFN storage
604 198 : IF (wf_history%store_wf_kp) THEN
605 104 : wf_history%store_wf = .FALSE.
606 104 : wf_history%store_overlap = .FALSE.
607 : ! store_wf_kp and store_overlap_kp remain TRUE
608 : ELSE
609 : ! Linear methods (except LINEAR_P) are still blocked
610 94 : IF (wf_history%store_wf .OR. wf_history%store_overlap) THEN
611 0 : CPABORT("Linear WFN-based extrapolation methods not implemented for k-points.")
612 : END IF
613 : END IF
614 198 : IF (wf_history%store_frozen_density) THEN
615 0 : CPABORT("Frozen density initialization method not possible for kpoints.")
616 : END IF
617 :
618 198 : END SUBROUTINE wfi_create_for_kp
619 :
620 : ! **************************************************************************************************
621 : !> \brief returns a string describing the interpolation method
622 : !> \param method_nr ...
623 : !> \return ...
624 : !> \par History
625 : !> 02.2003 created [fawzi]
626 : !> \author fawzi
627 : ! **************************************************************************************************
628 10815 : FUNCTION wfi_get_method_label(method_nr) RESULT(res)
629 : INTEGER, INTENT(in) :: method_nr
630 : CHARACTER(len=30) :: res
631 :
632 10815 : res = "unknown"
633 11053 : SELECT CASE (method_nr)
634 : CASE (wfi_use_prev_p_method_nr)
635 238 : res = "previous_p"
636 : CASE (wfi_use_prev_wf_method_nr)
637 213 : res = "previous_wf"
638 : CASE (wfi_use_prev_rho_r_method_nr)
639 4 : res = "previous_rho_r"
640 : CASE (wfi_use_guess_method_nr)
641 3442 : res = "initial_guess"
642 : CASE (wfi_linear_wf_method_nr)
643 2 : res = "mo linear"
644 : CASE (wfi_linear_p_method_nr)
645 3 : res = "P linear"
646 : CASE (wfi_linear_ps_method_nr)
647 6 : res = "PS linear"
648 : CASE (wfi_ps_method_nr)
649 188 : res = "PS Nth order"
650 : CASE (wfi_frozen_method_nr)
651 4 : res = "frozen density approximation"
652 : CASE (wfi_aspc_nr)
653 6691 : res = "ASPC"
654 : CASE (wfi_gext_proj_nr)
655 12 : res = "GEXT_PROJ"
656 : CASE (wfi_gext_proj_qtr_nr)
657 12 : res = "GEXT_PROJ_QTR"
658 : CASE default
659 : CALL cp_abort(__LOCATION__, &
660 : "Unknown interpolation method: "// &
661 10815 : TRIM(ADJUSTL(cp_to_string(method_nr))))
662 : END SELECT
663 10815 : END FUNCTION wfi_get_method_label
664 :
665 : ! **************************************************************************************************
666 : !> \brief calculates the new starting state for the scf for the next
667 : !> wf optimization
668 : !> \param wf_history the previous history needed to extrapolate
669 : !> \param qs_env the qs env with the latest result, and that will contain
670 : !> the new starting state
671 : !> \param dt the time at which to extrapolate (wrt. to the last snapshot)
672 : !> \param extrapolation_method_nr returns the extrapolation method used
673 : !> \param orthogonal_wf ...
674 : !> \par History
675 : !> 02.2003 created [fawzi]
676 : !> 11.2003 Joost VandeVondele : Implemented Nth order PS extrapolation
677 : !> 04.2026 Michele Nottoli : Added GEXT_PROJ and GEXT_PROJ_QTR extrapolations
678 : !> \author fawzi
679 : ! **************************************************************************************************
680 21295 : SUBROUTINE wfi_extrapolate(wf_history, qs_env, dt, extrapolation_method_nr, &
681 : orthogonal_wf)
682 : TYPE(qs_wf_history_type), POINTER :: wf_history
683 : TYPE(qs_environment_type), POINTER :: qs_env
684 : REAL(KIND=dp), INTENT(IN) :: dt
685 : INTEGER, INTENT(OUT), OPTIONAL :: extrapolation_method_nr
686 : LOGICAL, INTENT(OUT), OPTIONAL :: orthogonal_wf
687 :
688 : CHARACTER(len=*), PARAMETER :: routineN = 'wfi_extrapolate'
689 :
690 : INTEGER :: actual_extrapolation_method_nr, handle, &
691 : i, img, io_unit, ispin, k, n, nmo, &
692 : nvec, print_level
693 : LOGICAL :: do_kpoints, my_orthogonal_wf, use_overlap
694 : REAL(KIND=dp) :: alpha, t0, t1, t2
695 21295 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coeffs
696 21295 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools
697 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct, matrix_struct_new
698 : TYPE(cp_fm_type) :: csc, fm_tmp
699 : TYPE(cp_fm_type), POINTER :: mo_coeff
700 : TYPE(cp_logger_type), POINTER :: logger
701 21295 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao, rho_frozen_ao
702 21295 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
703 21295 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
704 : TYPE(qs_rho_type), POINTER :: rho
705 : TYPE(qs_wf_snapshot_type), POINTER :: t0_state, t1_state
706 :
707 21295 : NULLIFY (mos, ao_mo_fm_pools, t0_state, t1_state, mo_coeff, &
708 21295 : rho, rho_ao, rho_frozen_ao)
709 :
710 21295 : use_overlap = wf_history%store_overlap
711 :
712 21295 : CALL timeset(routineN, handle)
713 21295 : logger => cp_get_default_logger()
714 21295 : print_level = logger%iter_info%print_level
715 : io_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%SCF%PRINT%PROGRAM_RUN_INFO", &
716 21295 : extension=".scfLog")
717 :
718 21295 : CPASSERT(ASSOCIATED(wf_history))
719 21295 : CPASSERT(wf_history%ref_count > 0)
720 21295 : CPASSERT(ASSOCIATED(qs_env))
721 21295 : CALL get_qs_env(qs_env, mos=mos, rho=rho, do_kpoints=do_kpoints)
722 21295 : CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
723 : ! chooses the method for this extrapolation
724 21295 : IF (wf_history%snapshot_count < 1) THEN
725 : actual_extrapolation_method_nr = wfi_use_guess_method_nr
726 : ELSE
727 14708 : actual_extrapolation_method_nr = wf_history%interpolation_method_nr
728 : END IF
729 :
730 8 : SELECT CASE (actual_extrapolation_method_nr)
731 : CASE (wfi_linear_wf_method_nr)
732 8 : IF (wf_history%snapshot_count < 2) THEN
733 4 : actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
734 : END IF
735 : CASE (wfi_linear_p_method_nr)
736 12 : IF (wf_history%snapshot_count < 2) THEN
737 6 : actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
738 : END IF
739 : CASE (wfi_linear_ps_method_nr)
740 14708 : IF (wf_history%snapshot_count < 2) THEN
741 6 : actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
742 : END IF
743 : END SELECT
744 :
745 21295 : IF (PRESENT(extrapolation_method_nr)) &
746 21295 : extrapolation_method_nr = actual_extrapolation_method_nr
747 21295 : my_orthogonal_wf = .FALSE.
748 :
749 8 : SELECT CASE (actual_extrapolation_method_nr)
750 : CASE (wfi_frozen_method_nr)
751 8 : CPASSERT(.NOT. do_kpoints)
752 8 : t0_state => wfi_get_snapshot(wf_history, wf_index=1)
753 8 : CPASSERT(ASSOCIATED(t0_state%rho_frozen))
754 :
755 8 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
756 8 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
757 :
758 8 : CALL qs_rho_get(t0_state%rho_frozen, rho_ao=rho_frozen_ao)
759 8 : CALL qs_rho_get(rho, rho_ao=rho_ao)
760 16 : DO ispin = 1, SIZE(rho_frozen_ao)
761 : CALL dbcsr_copy(rho_ao(ispin)%matrix, &
762 : rho_frozen_ao(ispin)%matrix, &
763 16 : keep_sparsity=.TRUE.)
764 : END DO
765 : !FM updating rho_ao directly with t0_state%rho_ao would have the
766 : !FM wrong matrix structure
767 8 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
768 8 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
769 :
770 8 : my_orthogonal_wf = .FALSE.
771 : CASE (wfi_use_prev_rho_r_method_nr)
772 8 : t0_state => wfi_get_snapshot(wf_history, wf_index=1)
773 8 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
774 8 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
775 8 : IF (do_kpoints) THEN
776 0 : CPASSERT(ASSOCIATED(t0_state%rho_ao_kp))
777 0 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
778 0 : DO ispin = 1, SIZE(t0_state%rho_ao_kp, 1)
779 0 : DO img = 1, SIZE(t0_state%rho_ao_kp, 2)
780 0 : IF (img > SIZE(rho_ao_kp, 2)) THEN
781 0 : CPWARN("Change in cell neighborlist: might affect quality of initial guess")
782 : ELSE
783 : CALL dbcsr_copy(rho_ao_kp(ispin, img)%matrix, &
784 : t0_state%rho_ao_kp(ispin, img)%matrix, &
785 0 : keep_sparsity=.TRUE.)
786 : END IF
787 : END DO
788 : END DO
789 : ELSE
790 8 : CPASSERT(ASSOCIATED(t0_state%rho_ao))
791 8 : CALL qs_rho_get(rho, rho_ao=rho_ao)
792 16 : DO ispin = 1, SIZE(t0_state%rho_ao)
793 : CALL dbcsr_copy(rho_ao(ispin)%matrix, &
794 : t0_state%rho_ao(ispin)%matrix, &
795 16 : keep_sparsity=.TRUE.)
796 : END DO
797 : END IF
798 : ! Why is rho_g valid at this point ?
799 8 : CALL qs_rho_set(rho, rho_g_valid=.TRUE.)
800 :
801 : ! does nothing
802 : CASE (wfi_use_prev_wf_method_nr)
803 426 : my_orthogonal_wf = .TRUE.
804 426 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
805 426 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
806 :
807 426 : IF (do_kpoints) THEN
808 6 : CALL wfi_use_prev_wf_kp(qs_env, io_unit, print_level)
809 : ELSE
810 420 : CALL qs_rho_get(rho, rho_ao=rho_ao)
811 888 : DO ispin = 1, SIZE(mos)
812 468 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
813 468 : CALL reorthogonalize_vectors(qs_env, v_matrix=mo_coeff, n_col=nmo)
814 1356 : CALL calculate_density_matrix(mo_set=mos(ispin), density_matrix=rho_ao(ispin)%matrix)
815 : END DO
816 420 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
817 420 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
818 : END IF
819 :
820 : CASE (wfi_use_prev_p_method_nr)
821 476 : t0_state => wfi_get_snapshot(wf_history, wf_index=1)
822 476 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
823 476 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
824 476 : IF (do_kpoints) THEN
825 218 : CPASSERT(ASSOCIATED(t0_state%rho_ao_kp))
826 218 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
827 524 : DO ispin = 1, SIZE(t0_state%rho_ao_kp, 1)
828 31248 : DO img = 1, SIZE(t0_state%rho_ao_kp, 2)
829 31030 : IF (img > SIZE(rho_ao_kp, 2)) THEN
830 18 : CPWARN("Change in cell neighborlist: might affect quality of initial guess")
831 : ELSE
832 : CALL dbcsr_copy(rho_ao_kp(ispin, img)%matrix, &
833 : t0_state%rho_ao_kp(ispin, img)%matrix, &
834 30706 : keep_sparsity=.TRUE.)
835 : END IF
836 : END DO
837 : END DO
838 : ELSE
839 258 : CPASSERT(ASSOCIATED(t0_state%rho_ao))
840 258 : CALL qs_rho_get(rho, rho_ao=rho_ao)
841 646 : DO ispin = 1, SIZE(t0_state%rho_ao)
842 : CALL dbcsr_copy(rho_ao(ispin)%matrix, &
843 : t0_state%rho_ao(ispin)%matrix, &
844 646 : keep_sparsity=.TRUE.)
845 : END DO
846 : END IF
847 : !FM updating rho_ao directly with t0_state%rho_ao would have the
848 : !FM wrong matrix structure
849 476 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
850 476 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
851 :
852 : CASE (wfi_use_guess_method_nr)
853 : !FM more clean to do it here, but it
854 : !FM might need to read a file (restart) and thus globenv
855 : !FM I do not want globenv here, thus done by the caller
856 : !FM (btw. it also needs the eigensolver, and unless you relocate it
857 : !FM gives circular dependencies)
858 6829 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
859 6829 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
860 : CASE (wfi_linear_wf_method_nr)
861 4 : CPASSERT(.NOT. do_kpoints)
862 4 : t0_state => wfi_get_snapshot(wf_history, wf_index=2)
863 4 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
864 4 : CPASSERT(ASSOCIATED(t0_state))
865 4 : CPASSERT(ASSOCIATED(t1_state))
866 4 : CPASSERT(ASSOCIATED(t0_state%wf))
867 4 : CPASSERT(ASSOCIATED(t1_state%wf))
868 4 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
869 4 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
870 :
871 4 : my_orthogonal_wf = .TRUE.
872 4 : t0 = 0.0_dp
873 4 : t1 = t1_state%dt
874 4 : t2 = t1 + dt
875 4 : CALL qs_rho_get(rho, rho_ao=rho_ao)
876 8 : DO ispin = 1, SIZE(mos)
877 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, &
878 4 : nmo=nmo)
879 : CALL cp_fm_scale_and_add(alpha=0.0_dp, &
880 : matrix_a=mo_coeff, &
881 : matrix_b=t1_state%wf(ispin), &
882 4 : beta=(t2 - t0)/(t1 - t0))
883 : ! this copy should be unnecessary
884 : CALL cp_fm_scale_and_add(alpha=1.0_dp, &
885 : matrix_a=mo_coeff, &
886 4 : beta=(t1 - t2)/(t1 - t0), matrix_b=t0_state%wf(ispin))
887 : CALL reorthogonalize_vectors(qs_env, &
888 : v_matrix=mo_coeff, &
889 4 : n_col=nmo)
890 : CALL calculate_density_matrix(mo_set=mos(ispin), &
891 12 : density_matrix=rho_ao(ispin)%matrix)
892 : END DO
893 4 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
894 :
895 : CALL qs_ks_did_change(qs_env%ks_env, &
896 4 : rho_changed=.TRUE.)
897 : CASE (wfi_linear_p_method_nr)
898 6 : t0_state => wfi_get_snapshot(wf_history, wf_index=2)
899 6 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
900 6 : CPASSERT(ASSOCIATED(t0_state))
901 6 : CPASSERT(ASSOCIATED(t1_state))
902 6 : IF (do_kpoints) THEN
903 2 : CPASSERT(ASSOCIATED(t0_state%rho_ao_kp))
904 2 : CPASSERT(ASSOCIATED(t1_state%rho_ao_kp))
905 : ELSE
906 4 : CPASSERT(ASSOCIATED(t0_state%rho_ao))
907 4 : CPASSERT(ASSOCIATED(t1_state%rho_ao))
908 : END IF
909 6 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
910 6 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
911 :
912 6 : t0 = 0.0_dp
913 6 : t1 = t1_state%dt
914 6 : t2 = t1 + dt
915 6 : IF (do_kpoints) THEN
916 2 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
917 4 : DO ispin = 1, SIZE(rho_ao_kp, 1)
918 528 : DO img = 1, SIZE(rho_ao_kp, 2)
919 524 : IF (img > SIZE(t0_state%rho_ao_kp, 2) .OR. &
920 2 : img > SIZE(t1_state%rho_ao_kp, 2)) THEN
921 22 : CPWARN("Change in cell neighborlist: might affect quality of initial guess")
922 : ELSE
923 : CALL dbcsr_add(rho_ao_kp(ispin, img)%matrix, t1_state%rho_ao_kp(ispin, img)%matrix, &
924 502 : alpha_scalar=0.0_dp, beta_scalar=(t2 - t0)/(t1 - t0)) ! this copy should be unnecessary
925 : CALL dbcsr_add(rho_ao_kp(ispin, img)%matrix, t0_state%rho_ao_kp(ispin, img)%matrix, &
926 502 : alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0))
927 : END IF
928 : END DO
929 : END DO
930 : ELSE
931 4 : CALL qs_rho_get(rho, rho_ao=rho_ao)
932 8 : DO ispin = 1, SIZE(rho_ao)
933 : CALL dbcsr_add(rho_ao(ispin)%matrix, t1_state%rho_ao(ispin)%matrix, &
934 4 : alpha_scalar=0.0_dp, beta_scalar=(t2 - t0)/(t1 - t0)) ! this copy should be unnecessary
935 : CALL dbcsr_add(rho_ao(ispin)%matrix, t0_state%rho_ao(ispin)%matrix, &
936 8 : alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0))
937 : END DO
938 : END IF
939 6 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
940 6 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
941 : CASE (wfi_linear_ps_method_nr)
942 : ! wf not calculated, extract with PSC renormalized?
943 : ! use wf_linear?
944 12 : CPASSERT(.NOT. do_kpoints)
945 12 : t0_state => wfi_get_snapshot(wf_history, wf_index=2)
946 12 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
947 12 : CPASSERT(ASSOCIATED(t0_state))
948 12 : CPASSERT(ASSOCIATED(t1_state))
949 12 : CPASSERT(ASSOCIATED(t0_state%wf))
950 12 : CPASSERT(ASSOCIATED(t1_state%wf))
951 12 : IF (wf_history%store_overlap) THEN
952 4 : CPASSERT(ASSOCIATED(t0_state%overlap))
953 4 : CPASSERT(ASSOCIATED(t1_state%overlap))
954 : END IF
955 12 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
956 12 : IF (nvec >= wf_history%memory_depth) THEN
957 12 : IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. (qs_env%scf_control%eps_scf_hist /= 0)) THEN
958 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
959 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
960 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
961 12 : ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
962 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
963 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
964 12 : ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
965 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
966 : END IF
967 : END IF
968 :
969 12 : my_orthogonal_wf = .TRUE.
970 : ! use PS_2=2 PS_1-PS_0
971 : ! C_2 comes from using PS_2 as a projector acting on C_1
972 12 : CALL qs_rho_get(rho, rho_ao=rho_ao)
973 24 : DO ispin = 1, SIZE(mos)
974 12 : NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
975 12 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
976 : CALL cp_fm_get_info(mo_coeff, nrow_global=n, ncol_global=k, &
977 12 : matrix_struct=matrix_struct)
978 : CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, &
979 12 : nrow_global=k, ncol_global=k)
980 12 : CALL cp_fm_create(csc, matrix_struct_new)
981 12 : CALL cp_fm_struct_release(matrix_struct_new)
982 :
983 12 : IF (use_overlap) THEN
984 4 : CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), mo_coeff, k)
985 4 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), mo_coeff, 0.0_dp, csc)
986 : ELSE
987 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
988 8 : t1_state%wf(ispin), 0.0_dp, csc)
989 : END IF
990 12 : CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, mo_coeff)
991 12 : CALL cp_fm_release(csc)
992 12 : CALL cp_fm_scale_and_add(-1.0_dp, mo_coeff, 2.0_dp, t1_state%wf(ispin))
993 : CALL reorthogonalize_vectors(qs_env, &
994 : v_matrix=mo_coeff, &
995 12 : n_col=k)
996 : CALL calculate_density_matrix(mo_set=mos(ispin), &
997 48 : density_matrix=rho_ao(ispin)%matrix)
998 : END DO
999 12 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1000 12 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
1001 :
1002 : CASE (wfi_ps_method_nr)
1003 : ! figure out the actual number of vectors to use in the extrapolation:
1004 376 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
1005 376 : CPASSERT(nvec > 0)
1006 376 : IF (nvec >= wf_history%memory_depth) THEN
1007 178 : IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. (qs_env%scf_control%eps_scf_hist /= 0)) THEN
1008 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1009 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1010 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1011 178 : ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
1012 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1013 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1014 178 : ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
1015 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1016 : END IF
1017 : END IF
1018 :
1019 376 : IF (do_kpoints) THEN
1020 4 : CALL wfi_extrapolate_ps_aspc_kp(wf_history, qs_env, nvec, io_unit, print_level)
1021 4 : my_orthogonal_wf = .TRUE.
1022 : ELSE
1023 372 : my_orthogonal_wf = .TRUE.
1024 822 : DO ispin = 1, SIZE(mos)
1025 450 : NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
1026 450 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1027 : CALL cp_fm_get_info(mo_coeff, nrow_global=n, ncol_global=k, &
1028 450 : matrix_struct=matrix_struct)
1029 450 : CALL cp_fm_create(fm_tmp, matrix_struct)
1030 : CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, &
1031 450 : nrow_global=k, ncol_global=k)
1032 450 : CALL cp_fm_create(csc, matrix_struct_new)
1033 450 : CALL cp_fm_struct_release(matrix_struct_new)
1034 : ! first the most recent
1035 450 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
1036 450 : CALL cp_fm_to_fm(t1_state%wf(ispin), mo_coeff)
1037 450 : alpha = nvec
1038 450 : CALL cp_fm_scale(alpha, mo_coeff)
1039 450 : CALL qs_rho_get(rho, rho_ao=rho_ao)
1040 962 : DO i = 2, nvec
1041 512 : t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1042 512 : IF (use_overlap) THEN
1043 474 : CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
1044 474 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
1045 : ELSE
1046 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
1047 38 : t1_state%wf(ispin), 0.0_dp, csc)
1048 : END IF
1049 512 : CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
1050 512 : alpha = -1.0_dp*alpha*REAL(nvec - i + 1, dp)/REAL(i, dp)
1051 962 : CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp)
1052 : END DO
1053 :
1054 450 : CALL cp_fm_release(csc)
1055 450 : CALL cp_fm_release(fm_tmp)
1056 : CALL reorthogonalize_vectors(qs_env, &
1057 : v_matrix=mo_coeff, &
1058 450 : n_col=k)
1059 : CALL calculate_density_matrix(mo_set=mos(ispin), &
1060 1722 : density_matrix=rho_ao(ispin)%matrix)
1061 : END DO
1062 372 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1063 372 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
1064 : END IF
1065 :
1066 : CASE (wfi_aspc_nr)
1067 13102 : CALL cite_reference(Kolafa2004)
1068 13102 : CALL cite_reference(Kuhne2007)
1069 : ! figure out the actual number of vectors to use in the extrapolation:
1070 13102 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
1071 13102 : CPASSERT(nvec > 0)
1072 13102 : IF (nvec >= wf_history%memory_depth) THEN
1073 8366 : IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. &
1074 : (qs_env%scf_control%eps_scf_hist /= 0)) THEN
1075 18 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1076 18 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1077 18 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1078 8348 : ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
1079 62 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1080 62 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1081 8286 : ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
1082 8 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1083 : END IF
1084 : END IF
1085 :
1086 13102 : IF (do_kpoints) THEN
1087 358 : CALL wfi_extrapolate_ps_aspc_kp(wf_history, qs_env, nvec, io_unit, print_level)
1088 358 : my_orthogonal_wf = .TRUE.
1089 : ELSE
1090 12744 : my_orthogonal_wf = .TRUE.
1091 12744 : CALL qs_rho_get(rho, rho_ao=rho_ao)
1092 26071 : DO ispin = 1, SIZE(mos)
1093 13327 : NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
1094 13327 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1095 : CALL cp_fm_get_info(mo_coeff, &
1096 : nrow_global=n, &
1097 : ncol_global=k, &
1098 13327 : matrix_struct=matrix_struct)
1099 13327 : CALL cp_fm_create(fm_tmp, matrix_struct, set_zero=.TRUE.)
1100 : CALL cp_fm_struct_create(matrix_struct_new, &
1101 : template_fmstruct=matrix_struct, &
1102 : nrow_global=k, &
1103 13327 : ncol_global=k)
1104 13327 : CALL cp_fm_create(csc, matrix_struct_new, set_zero=.TRUE.)
1105 13327 : CALL cp_fm_struct_release(matrix_struct_new)
1106 : ! first the most recent
1107 : t1_state => wfi_get_snapshot(wf_history, &
1108 13327 : wf_index=1)
1109 13327 : CALL cp_fm_to_fm(t1_state%wf(ispin), mo_coeff)
1110 13327 : alpha = REAL(4*nvec - 2, KIND=dp)/REAL(nvec + 1, KIND=dp)
1111 13327 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1112 : WRITE (UNIT=io_unit, FMT="(/,T2,A,/,/,T3,A,I0,/,/,T3,A2,I0,A4,F10.6)") &
1113 2505 : "Parameters for the always stable predictor-corrector (ASPC) method:", &
1114 2505 : "ASPC order: ", MAX(nvec - 2, 0), &
1115 5010 : "B(", 1, ") = ", alpha
1116 : END IF
1117 13327 : CALL cp_fm_scale(alpha, mo_coeff)
1118 :
1119 51237 : DO i = 2, nvec
1120 37910 : t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1121 37910 : IF (use_overlap) THEN
1122 26370 : CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
1123 26370 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
1124 : ELSE
1125 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
1126 11540 : t1_state%wf(ispin), 0.0_dp, csc)
1127 : END IF
1128 37910 : CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
1129 : alpha = (-1.0_dp)**(i + 1)*REAL(i, KIND=dp)* &
1130 37910 : binomial(2*nvec, nvec - i)/binomial(2*nvec - 2, nvec - 1)
1131 37910 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1132 : WRITE (UNIT=io_unit, FMT="(T3,A2,I0,A4,F10.6)") &
1133 7212 : "B(", i, ") = ", alpha
1134 : END IF
1135 51237 : CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp)
1136 : END DO
1137 13327 : CALL cp_fm_release(csc)
1138 13327 : CALL cp_fm_release(fm_tmp)
1139 : CALL reorthogonalize_vectors(qs_env, &
1140 : v_matrix=mo_coeff, &
1141 13327 : n_col=k)
1142 : CALL calculate_density_matrix(mo_set=mos(ispin), &
1143 39398 : density_matrix=rho_ao(ispin)%matrix)
1144 : END DO
1145 12744 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1146 12744 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
1147 : END IF ! do_kpoints
1148 :
1149 : CASE (wfi_gext_proj_nr)
1150 24 : CPASSERT(.NOT. do_kpoints)
1151 :
1152 : ! figure out the actual number of vectors to use in the extrapolation:
1153 24 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
1154 24 : IF (nvec >= wf_history%memory_depth) THEN
1155 8 : IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. &
1156 : (qs_env%scf_control%eps_scf_hist /= 0)) THEN
1157 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1158 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1159 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1160 8 : ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
1161 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1162 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1163 8 : ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
1164 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1165 : END IF
1166 : END IF
1167 24 : CPASSERT(nvec > 0)
1168 :
1169 : ! get the coefficients for the fitting
1170 72 : ALLOCATE (coeffs(nvec))
1171 24 : NULLIFY (matrix_s)
1172 24 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
1173 : CALL diff_fitting(wf_history, matrix_s(1)%matrix, coeffs, nvec, &
1174 24 : 1e-4_dp, io_unit, print_level)
1175 :
1176 24 : my_orthogonal_wf = .TRUE.
1177 24 : CALL qs_rho_get(rho, rho_ao=rho_ao)
1178 48 : DO ispin = 1, SIZE(mos)
1179 24 : NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
1180 24 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1181 : CALL cp_fm_get_info(mo_coeff, &
1182 : nrow_global=n, &
1183 : ncol_global=k, &
1184 24 : matrix_struct=matrix_struct)
1185 24 : CALL cp_fm_create(fm_tmp, matrix_struct)
1186 : CALL cp_fm_struct_create(matrix_struct_new, &
1187 : template_fmstruct=matrix_struct, &
1188 : nrow_global=k, &
1189 24 : ncol_global=k)
1190 24 : CALL cp_fm_create(csc, matrix_struct_new)
1191 24 : CALL cp_fm_struct_release(matrix_struct_new)
1192 :
1193 24 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
1194 :
1195 : ! do the linear combination of previous PSs
1196 24 : CALL cp_fm_set_all(mo_coeff, 0.0_dp)
1197 104 : DO i = 1, nvec
1198 80 : t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1199 80 : CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
1200 80 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
1201 80 : CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
1202 104 : CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, coeffs(i), fm_tmp)
1203 : END DO
1204 24 : CALL cp_fm_release(csc)
1205 24 : CALL cp_fm_release(fm_tmp)
1206 : CALL reorthogonalize_vectors(qs_env, &
1207 : v_matrix=mo_coeff, &
1208 24 : n_col=k)
1209 : CALL calculate_density_matrix(mo_set=mos(ispin), &
1210 96 : density_matrix=rho_ao(ispin)%matrix)
1211 : END DO
1212 24 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1213 24 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
1214 :
1215 24 : DEALLOCATE (coeffs)
1216 :
1217 : CASE (wfi_gext_proj_qtr_nr)
1218 24 : CPASSERT(.NOT. do_kpoints)
1219 :
1220 : ! figure out the actual number of vectors to use in the extrapolation:
1221 24 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
1222 24 : IF (nvec >= wf_history%memory_depth) THEN
1223 8 : IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. &
1224 : (qs_env%scf_control%eps_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%eps_scf = qs_env%scf_control%eps_scf_hist
1227 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1228 8 : ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
1229 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1230 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1231 8 : ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
1232 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1233 : END IF
1234 : END IF
1235 24 : CPASSERT(nvec > 0)
1236 :
1237 : ! get the coefficients for the fitting
1238 72 : ALLOCATE (coeffs(nvec))
1239 24 : NULLIFY (matrix_s)
1240 24 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
1241 : CALL tr_fitting(wf_history, matrix_s(1)%matrix, coeffs, nvec, &
1242 24 : 1e-4_dp, io_unit, print_level)
1243 :
1244 24 : my_orthogonal_wf = .TRUE.
1245 24 : CALL qs_rho_get(rho, rho_ao=rho_ao)
1246 48 : DO ispin = 1, SIZE(mos)
1247 24 : NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
1248 24 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1249 : CALL cp_fm_get_info(mo_coeff, &
1250 : nrow_global=n, &
1251 : ncol_global=k, &
1252 24 : matrix_struct=matrix_struct)
1253 24 : CALL cp_fm_create(fm_tmp, matrix_struct)
1254 : CALL cp_fm_struct_create(matrix_struct_new, &
1255 : template_fmstruct=matrix_struct, &
1256 : nrow_global=k, &
1257 24 : ncol_global=k)
1258 24 : CALL cp_fm_create(csc, matrix_struct_new)
1259 24 : CALL cp_fm_struct_release(matrix_struct_new)
1260 :
1261 24 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
1262 :
1263 : ! do the linear combination of previous PSs
1264 24 : CALL cp_fm_set_all(mo_coeff, 0.0_dp)
1265 104 : DO i = 1, nvec
1266 80 : t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1267 80 : CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
1268 80 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
1269 80 : CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
1270 104 : CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, coeffs(i), fm_tmp)
1271 : END DO
1272 24 : CALL cp_fm_release(csc)
1273 24 : CALL cp_fm_release(fm_tmp)
1274 : CALL reorthogonalize_vectors(qs_env, &
1275 : v_matrix=mo_coeff, &
1276 24 : n_col=k)
1277 : CALL calculate_density_matrix(mo_set=mos(ispin), &
1278 96 : density_matrix=rho_ao(ispin)%matrix)
1279 : END DO
1280 24 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1281 24 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
1282 :
1283 24 : DEALLOCATE (coeffs)
1284 :
1285 : CASE default
1286 : CALL cp_abort(__LOCATION__, &
1287 : "Unknown interpolation method: "// &
1288 21295 : TRIM(ADJUSTL(cp_to_string(wf_history%interpolation_method_nr))))
1289 : END SELECT
1290 21295 : IF (PRESENT(orthogonal_wf)) orthogonal_wf = my_orthogonal_wf
1291 : CALL cp_print_key_finished_output(io_unit, logger, qs_env%input, &
1292 21295 : "DFT%SCF%PRINT%PROGRAM_RUN_INFO")
1293 21295 : CALL timestop(handle)
1294 21295 : END SUBROUTINE wfi_extrapolate
1295 :
1296 : ! **************************************************************************************************
1297 : !> \brief Reorthogonalizes the wavefunctions from the previous step for k-points
1298 : !> using the current S(k) metric and rebuilds the density matrix.
1299 : !> \param qs_env The QS environment
1300 : !> \param io_unit output unit
1301 : !> \param print_level print level
1302 : ! **************************************************************************************************
1303 368 : SUBROUTINE wfi_use_prev_wf_kp(qs_env, io_unit, print_level)
1304 : TYPE(qs_environment_type), POINTER :: qs_env
1305 : INTEGER, INTENT(IN) :: io_unit, print_level
1306 :
1307 : CHARACTER(len=*), PARAMETER :: routineN = 'wfi_use_prev_wf_kp'
1308 :
1309 368 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: col_scaling
1310 : INTEGER :: chol_info, handle, igroup, ik, ikp, &
1311 : indx, ispin, j, kplocal, nao, nkp, &
1312 : nkp_groups, nmo, nspin
1313 : INTEGER, DIMENSION(2) :: kp_range
1314 368 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
1315 368 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1316 : LOGICAL :: my_kpgrp, use_real_wfn
1317 : REAL(KIND=dp) :: eval_thresh
1318 368 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1319 368 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1320 368 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
1321 : TYPE(cp_cfm_type) :: cfm_evecs, cfm_mhalf, cfm_nao_nmo_work, &
1322 : cmos_new, csc_cfm
1323 368 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: csmat_cur
1324 368 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools_kp
1325 : TYPE(cp_fm_struct_type), POINTER :: ao_ao_struct, nmo_nmo_struct
1326 : TYPE(cp_fm_type) :: fmdummy, fmlocal
1327 : TYPE(cp_fm_type), POINTER :: imos, mo_coeff, rmos
1328 368 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, rho_ao_kp
1329 : TYPE(dbcsr_type), POINTER :: cmatrix_db, rmatrix, tmpmat
1330 : TYPE(dft_control_type), POINTER :: dft_control
1331 : TYPE(kpoint_env_type), POINTER :: kp
1332 : TYPE(kpoint_type), POINTER :: kpoints
1333 : TYPE(mp_para_env_type), POINTER :: para_env
1334 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1335 368 : POINTER :: sab_nl
1336 : TYPE(qs_matrix_pools_type), POINTER :: mpools_kp
1337 : TYPE(qs_rho_type), POINTER :: rho
1338 : TYPE(qs_scf_env_type), POINTER :: scf_env
1339 : TYPE(scf_control_type), POINTER :: scf_control
1340 :
1341 368 : CALL timeset(routineN, handle)
1342 :
1343 368 : NULLIFY (kpoints, dft_control, matrix_s_kp, scf_env, scf_control, rho, sab_nl, kp, mo_coeff, rmos, imos)
1344 :
1345 : CALL get_qs_env(qs_env, kpoints=kpoints, dft_control=dft_control, &
1346 : matrix_s_kp=matrix_s_kp, scf_env=scf_env, &
1347 368 : scf_control=scf_control, rho=rho)
1348 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn, &
1349 : kp_range=kp_range, nkp_groups=nkp_groups, kp_dist=kp_dist, &
1350 368 : sab_nl=sab_nl, cell_to_index=cell_to_index)
1351 368 : kplocal = kp_range(2) - kp_range(1) + 1
1352 :
1353 368 : IF (use_real_wfn) THEN
1354 0 : CALL timestop(handle)
1355 0 : RETURN
1356 : END IF
1357 :
1358 368 : kp => kpoints%kp_env(1)%kpoint_env
1359 368 : nspin = SIZE(kp%mos, 2)
1360 368 : CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo, mo_coeff=mo_coeff)
1361 :
1362 368 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1363 : WRITE (UNIT=io_unit, FMT="(/,T2,A)") &
1364 0 : "Using previous wavefunctions as initial guess for k-points (with reorthogonalization)"
1365 : END IF
1366 :
1367 : ! Allocate dbcsr work matrices
1368 368 : ALLOCATE (rmatrix, cmatrix_db, tmpmat)
1369 368 : CALL dbcsr_create(rmatrix, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
1370 368 : CALL dbcsr_create(cmatrix_db, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
1371 368 : CALL dbcsr_create(tmpmat, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1372 368 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
1373 368 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix_db, sab_nl)
1374 :
1375 368 : CALL get_kpoint_info(kpoints, mpools=mpools_kp)
1376 368 : CALL mpools_get(mpools_kp, ao_ao_fm_pools=ao_ao_fm_pools_kp)
1377 368 : CALL fm_pool_create_fm(ao_ao_fm_pools_kp(1)%pool, fmlocal)
1378 368 : CALL cp_fm_get_info(fmlocal, matrix_struct=ao_ao_struct)
1379 :
1380 368 : CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct)
1381 368 : CALL cp_cfm_create(cfm_nao_nmo_work, mo_coeff%matrix_struct)
1382 :
1383 368 : NULLIFY (nmo_nmo_struct)
1384 : CALL cp_fm_struct_create(nmo_nmo_struct, template_fmstruct=mo_coeff%matrix_struct, &
1385 368 : nrow_global=nmo, ncol_global=nmo)
1386 368 : CALL cp_cfm_create(csc_cfm, nmo_nmo_struct)
1387 368 : CALL cp_fm_struct_release(nmo_nmo_struct)
1388 :
1389 368 : para_env => kpoints%blacs_env_all%para_env
1390 7456 : ALLOCATE (info(kplocal*nkp_groups, 2))
1391 :
1392 2104 : ALLOCATE (csmat_cur(kplocal))
1393 1368 : DO ikp = 1, kplocal
1394 1368 : CALL cp_cfm_create(csmat_cur(ikp), ao_ao_struct)
1395 : END DO
1396 :
1397 : ! Phase A: Fourier Transform S(R) -> S(k)
1398 : indx = 0
1399 1368 : DO ikp = 1, kplocal
1400 3072 : DO igroup = 1, nkp_groups
1401 1704 : ik = kp_dist(1, igroup) + ikp - 1
1402 1704 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
1403 1704 : indx = indx + 1
1404 :
1405 1704 : CALL dbcsr_set(rmatrix, 0.0_dp)
1406 1704 : CALL dbcsr_set(cmatrix_db, 0.0_dp)
1407 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix_db, rsmat=matrix_s_kp, &
1408 1704 : ispin=1, xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1409 1704 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1410 1704 : CALL copy_dbcsr_to_fm(tmpmat, scf_env%scf_work1(1))
1411 1704 : CALL dbcsr_desymmetrize(cmatrix_db, tmpmat)
1412 1704 : CALL copy_dbcsr_to_fm(tmpmat, scf_env%scf_work1(2))
1413 :
1414 2704 : IF (my_kpgrp) THEN
1415 1000 : CALL cp_fm_start_copy_general(scf_env%scf_work1(1), fmlocal, para_env, info(indx, 1))
1416 1000 : CALL cp_fm_start_copy_general(scf_env%scf_work1(2), fmlocal, para_env, info(indx, 2))
1417 : ELSE
1418 704 : CALL cp_fm_start_copy_general(scf_env%scf_work1(1), fmdummy, para_env, info(indx, 1))
1419 704 : CALL cp_fm_start_copy_general(scf_env%scf_work1(2), fmdummy, para_env, info(indx, 2))
1420 : END IF
1421 : END DO
1422 : END DO
1423 :
1424 : ! Finish Communication
1425 : indx = 0
1426 1368 : DO ikp = 1, kplocal
1427 3072 : DO igroup = 1, nkp_groups
1428 1704 : ik = kp_dist(1, igroup) + ikp - 1
1429 1704 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
1430 704 : indx = indx + 1
1431 1000 : IF (my_kpgrp) THEN
1432 1000 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
1433 1000 : CALL cp_cfm_scale_and_add_fm(z_zero, csmat_cur(ikp), z_one, fmlocal)
1434 1000 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
1435 1000 : CALL cp_cfm_scale_and_add_fm(z_one, csmat_cur(ikp), gaussi, fmlocal)
1436 : END IF
1437 : END DO
1438 : END DO
1439 :
1440 2072 : DO indx = 1, kplocal*nkp_groups
1441 1704 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
1442 2072 : CALL cp_fm_cleanup_copy_general(info(indx, 2))
1443 : END DO
1444 :
1445 : ! Phase B: Orthogonalize current MOs w.r.t. new S(k)
1446 1104 : ALLOCATE (eigenvalues(nmo))
1447 368 : eval_thresh = 1.0E-12_dp
1448 :
1449 1368 : DO ikp = 1, kplocal
1450 1000 : kp => kpoints%kp_env(ikp)%kpoint_env
1451 2560 : DO ispin = 1, nspin
1452 1192 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
1453 1192 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1454 1192 : CALL cp_fm_to_cfm(rmos, imos, cmos_new)
1455 :
1456 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nao, z_one, &
1457 1192 : csmat_cur(ikp), cmos_new, z_zero, cfm_nao_nmo_work)
1458 : CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, &
1459 1192 : cmos_new, cfm_nao_nmo_work, z_zero, csc_cfm)
1460 :
1461 1192 : CALL cp_cfm_cholesky_decompose(csc_cfm, info_out=chol_info)
1462 1192 : IF (chol_info == 0) THEN
1463 1192 : CALL cp_cfm_triangular_multiply(csc_cfm, cmos_new, side='R', invert_tr=.TRUE., uplo_tr='U')
1464 : ELSE
1465 0 : CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, cmos_new, cfm_nao_nmo_work, z_zero, csc_cfm)
1466 0 : CALL cp_cfm_create(cfm_evecs, csc_cfm%matrix_struct)
1467 0 : CALL cp_cfm_create(cfm_mhalf, csc_cfm%matrix_struct)
1468 0 : CALL cp_cfm_heevd(csc_cfm, cfm_evecs, eigenvalues)
1469 0 : CALL cp_cfm_to_cfm(cfm_evecs, cfm_mhalf)
1470 0 : ALLOCATE (col_scaling(nmo))
1471 0 : DO j = 1, nmo
1472 0 : IF (eigenvalues(j) > eval_thresh) THEN
1473 0 : col_scaling(j) = CMPLX(1.0_dp/SQRT(eigenvalues(j)), 0.0_dp, KIND=dp)
1474 : ELSE
1475 0 : col_scaling(j) = z_zero
1476 : END IF
1477 : END DO
1478 0 : CALL cp_cfm_column_scale(cfm_mhalf, col_scaling)
1479 0 : DEALLOCATE (col_scaling)
1480 0 : CALL cp_cfm_gemm('N', 'C', nmo, nmo, nmo, z_one, cfm_mhalf, cfm_evecs, z_zero, csc_cfm)
1481 0 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nmo, z_one, cmos_new, csc_cfm, z_zero, cfm_nao_nmo_work)
1482 0 : CALL cp_cfm_to_cfm(cfm_nao_nmo_work, cmos_new)
1483 0 : CALL cp_cfm_release(cfm_evecs)
1484 0 : CALL cp_cfm_release(cfm_mhalf)
1485 : END IF
1486 3384 : CALL cp_cfm_to_fm(cmos_new, rmos, imos)
1487 : END DO
1488 : END DO
1489 368 : DEALLOCATE (eigenvalues)
1490 :
1491 : ! Phase C: Rebuild Density Matrix P(R)
1492 368 : CALL kpoint_set_mo_occupation(kpoints, scf_control%smear)
1493 368 : CALL kpoint_density_matrices(kpoints)
1494 368 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1495 : CALL kpoint_density_transform(kpoints, rho_ao_kp, .FALSE., &
1496 368 : matrix_s_kp(1, 1)%matrix, sab_nl, scf_env%scf_work1)
1497 368 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1498 368 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
1499 :
1500 : ! Cleanup
1501 1368 : DO ikp = 1, kplocal
1502 1368 : CALL cp_cfm_release(csmat_cur(ikp))
1503 : END DO
1504 368 : DEALLOCATE (csmat_cur)
1505 3776 : DEALLOCATE (info)
1506 368 : CALL cp_cfm_release(cmos_new)
1507 368 : CALL cp_cfm_release(cfm_nao_nmo_work)
1508 368 : CALL cp_cfm_release(csc_cfm)
1509 368 : CALL fm_pool_give_back_fm(ao_ao_fm_pools_kp(1)%pool, fmlocal)
1510 368 : CALL dbcsr_deallocate_matrix(rmatrix)
1511 368 : CALL dbcsr_deallocate_matrix(cmatrix_db)
1512 368 : CALL dbcsr_deallocate_matrix(tmpmat)
1513 :
1514 368 : CALL timestop(handle)
1515 2208 : END SUBROUTINE wfi_use_prev_wf_kp
1516 :
1517 : ! **************************************************************************************************
1518 : !> \brief Performs PS/ASPC wavefunction extrapolation for k-point calculations.
1519 : !> Applies PS/ASPC coefficients to complex MO coefficients at each k-point,
1520 : !> with subspace alignment via historical overlap matrices.
1521 : !> Delegates final orthogonalization and density building to wfi_use_prev_wf_kp.
1522 : !> \param wf_history wavefunction history buffer
1523 : !> \param qs_env QS environment
1524 : !> \param nvec number of history snapshots to use
1525 : !> \param io_unit output unit for logging
1526 : !> \param print_level current print level
1527 : ! **************************************************************************************************
1528 724 : SUBROUTINE wfi_extrapolate_ps_aspc_kp(wf_history, qs_env, nvec, io_unit, print_level)
1529 : TYPE(qs_wf_history_type), POINTER :: wf_history
1530 : TYPE(qs_environment_type), POINTER :: qs_env
1531 : INTEGER, INTENT(IN) :: nvec, io_unit, print_level
1532 :
1533 : CHARACTER(len=*), PARAMETER :: routineN = 'wfi_extrapolate_ps_aspc_kp'
1534 :
1535 : INTEGER :: handle, i, ikp, ispin, kplocal, &
1536 : method_nr, nao, nmo, nspin
1537 : INTEGER, DIMENSION(2) :: kp_range
1538 : LOGICAL :: use_real_wfn
1539 : REAL(KIND=dp) :: alpha_coeff
1540 : TYPE(cp_cfm_type) :: cfm_nao_nmo_work, cmos_1, cmos_i, &
1541 : cmos_new, csc_cfm
1542 : TYPE(cp_fm_struct_type), POINTER :: nmo_nmo_struct
1543 : TYPE(cp_fm_type), POINTER :: imos, mo_coeff, rmos
1544 : TYPE(kpoint_env_type), POINTER :: kp
1545 : TYPE(kpoint_type), POINTER :: kpoints
1546 : TYPE(qs_wf_snapshot_type), POINTER :: t0_state, t1_state
1547 :
1548 362 : method_nr = wf_history%interpolation_method_nr
1549 :
1550 362 : CALL timeset(routineN, handle)
1551 362 : NULLIFY (kpoints, kp, mo_coeff, rmos, imos, t0_state, t1_state, nmo_nmo_struct)
1552 :
1553 362 : CALL get_qs_env(qs_env, kpoints=kpoints)
1554 362 : CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn, kp_range=kp_range)
1555 362 : kplocal = kp_range(2) - kp_range(1) + 1
1556 :
1557 362 : IF (use_real_wfn) THEN
1558 0 : IF (method_nr == wfi_aspc_nr) THEN
1559 : CALL cp_warn(__LOCATION__, "ASPC with k-points requires complex wavefunctions; "// &
1560 0 : "falling back to USE_PREV_WF.")
1561 : ELSE
1562 : CALL cp_warn(__LOCATION__, "PS with k-points requires complex wavefunctions; "// &
1563 0 : "falling back to USE_PREV_WF.")
1564 : END IF
1565 0 : CALL wfi_use_prev_wf_kp(qs_env, io_unit, print_level)
1566 0 : CALL timestop(handle)
1567 0 : RETURN
1568 : END IF
1569 :
1570 362 : kp => kpoints%kp_env(1)%kpoint_env
1571 362 : nspin = SIZE(kp%mos, 2)
1572 362 : CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo, mo_coeff=mo_coeff)
1573 :
1574 362 : IF (method_nr == wfi_aspc_nr) THEN
1575 358 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1576 : WRITE (UNIT=io_unit, FMT="(/,T2,A,/,T3,A,I0)") &
1577 7 : "Parameters for the always stable predictor-corrector (ASPC) method:", &
1578 14 : "ASPC order: ", MAX(nvec - 2, 0)
1579 : END IF
1580 : END IF
1581 :
1582 16 : IF (method_nr == wfi_aspc_nr) THEN
1583 358 : CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct, set_zero=.TRUE.)
1584 358 : CALL cp_cfm_create(cmos_1, mo_coeff%matrix_struct, set_zero=.TRUE.)
1585 358 : CALL cp_cfm_create(cmos_i, mo_coeff%matrix_struct, set_zero=.TRUE.)
1586 358 : CALL cp_cfm_create(cfm_nao_nmo_work, mo_coeff%matrix_struct, set_zero=.TRUE.)
1587 : ELSE
1588 4 : CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct)
1589 4 : CALL cp_cfm_create(cmos_1, mo_coeff%matrix_struct)
1590 4 : CALL cp_cfm_create(cmos_i, mo_coeff%matrix_struct)
1591 4 : CALL cp_cfm_create(cfm_nao_nmo_work, mo_coeff%matrix_struct)
1592 : END IF
1593 :
1594 : CALL cp_fm_struct_create(nmo_nmo_struct, template_fmstruct=mo_coeff%matrix_struct, &
1595 362 : nrow_global=nmo, ncol_global=nmo)
1596 362 : IF (method_nr == wfi_aspc_nr) THEN
1597 358 : CALL cp_cfm_create(csc_cfm, nmo_nmo_struct, set_zero=.TRUE.)
1598 : ELSE
1599 4 : CALL cp_cfm_create(csc_cfm, nmo_nmo_struct)
1600 : END IF
1601 362 : CALL cp_fm_struct_release(nmo_nmo_struct)
1602 :
1603 : ! Phase 1: Initialize C_new(k) = B(1) * C_1(k)
1604 362 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
1605 362 : IF (method_nr == wfi_aspc_nr) THEN
1606 358 : alpha_coeff = REAL(4*nvec - 2, KIND=dp)/REAL(nvec + 1, KIND=dp)
1607 358 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1608 7 : WRITE (UNIT=io_unit, FMT="(T3,A2,I0,A4,F10.6)") "B(", 1, ") = ", alpha_coeff
1609 : END IF
1610 : ELSE
1611 4 : alpha_coeff = nvec
1612 : END IF
1613 :
1614 1338 : DO ikp = 1, kplocal
1615 976 : kp => kpoints%kp_env(ikp)%kpoint_env
1616 2506 : DO ispin = 1, nspin
1617 1168 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
1618 1168 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1619 1168 : CALL cp_fm_to_fm(t1_state%wf_kp(ikp, 1, ispin), rmos)
1620 1168 : CALL cp_fm_to_fm(t1_state%wf_kp(ikp, 2, ispin), imos)
1621 1168 : CALL cp_fm_scale(alpha_coeff, rmos)
1622 2144 : CALL cp_fm_scale(alpha_coeff, imos)
1623 : END DO
1624 : END DO
1625 :
1626 : ! Phase 2: Accumulate historical snapshots C_new += B(i) * C_proj(k)
1627 1544 : DO i = 2, nvec
1628 1182 : t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1629 1182 : IF (method_nr == wfi_aspc_nr) THEN
1630 : alpha_coeff = (-1.0_dp)**(i + 1)*REAL(i, KIND=dp)* &
1631 1180 : binomial(2*nvec, nvec - i)/binomial(2*nvec - 2, nvec - 1)
1632 1180 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1633 5 : WRITE (UNIT=io_unit, FMT="(T3,A2,I0,A4,F10.6)") "B(", i, ") = ", alpha_coeff
1634 : END IF
1635 : ELSE
1636 2 : alpha_coeff = -1.0_dp*alpha_coeff*REAL(nvec - i + 1, dp)/REAL(i, dp)
1637 : END IF
1638 :
1639 3816 : DO ikp = 1, kplocal
1640 2272 : kp => kpoints%kp_env(ikp)%kpoint_env
1641 5822 : DO ispin = 1, nspin
1642 2368 : CALL cp_fm_to_cfm(t1_state%wf_kp(ikp, 1, ispin), t1_state%wf_kp(ikp, 2, ispin), cmos_1)
1643 2368 : CALL cp_fm_to_cfm(t0_state%wf_kp(ikp, 1, ispin), t0_state%wf_kp(ikp, 2, ispin), cmos_i)
1644 :
1645 : ! Subspace projection: C_proj = C_i * (C_i^dag S_i C_1)
1646 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nao, z_one, &
1647 2368 : t0_state%overlap_cfm_kp(ikp), cmos_1, z_zero, cfm_nao_nmo_work)
1648 : CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, &
1649 2368 : cmos_i, cfm_nao_nmo_work, z_zero, csc_cfm)
1650 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nmo, z_one, &
1651 2368 : cmos_i, csc_cfm, z_zero, cfm_nao_nmo_work)
1652 :
1653 2368 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
1654 2368 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1655 2368 : CALL cp_fm_to_cfm(rmos, imos, cmos_new)
1656 2368 : CALL cp_cfm_scale_and_add(z_one, cmos_new, CMPLX(alpha_coeff, 0.0_dp, KIND=dp), cfm_nao_nmo_work)
1657 4640 : CALL cp_cfm_to_fm(cmos_new, rmos, imos)
1658 : END DO
1659 : END DO
1660 : END DO
1661 :
1662 362 : CALL cp_cfm_release(cmos_new)
1663 362 : CALL cp_cfm_release(cmos_1)
1664 362 : CALL cp_cfm_release(cmos_i)
1665 362 : CALL cp_cfm_release(cfm_nao_nmo_work)
1666 362 : CALL cp_cfm_release(csc_cfm)
1667 :
1668 : ! Phase 3: Delegate Orthogonalization and Density Building
1669 : ! We pass io_unit = 0 to suppress the redundant "Using previous..." print
1670 : ! inside wfi_use_prev_wf_kp, keeping the ASPC log output perfectly clean.
1671 362 : CALL wfi_use_prev_wf_kp(qs_env, 0, print_level)
1672 :
1673 362 : CALL timestop(handle)
1674 :
1675 362 : END SUBROUTINE wfi_extrapolate_ps_aspc_kp
1676 :
1677 : ! **************************************************************************************************
1678 : !> \brief Decides if scf control variables has to changed due
1679 : !> to using a WF extrapolation.
1680 : !> \param qs_env The QS environment
1681 : !> \param nvec ...
1682 : !> \par History
1683 : !> 11.2006 created [TdK]
1684 : !> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch)
1685 : ! **************************************************************************************************
1686 7757 : ELEMENTAL SUBROUTINE wfi_set_history_variables(qs_env, nvec)
1687 : TYPE(qs_environment_type), INTENT(INOUT) :: qs_env
1688 : INTEGER, INTENT(IN) :: nvec
1689 :
1690 7757 : IF (nvec >= qs_env%wf_history%memory_depth) THEN
1691 1289 : IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. (qs_env%scf_control%eps_scf_hist /= 0)) THEN
1692 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1693 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1694 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1695 1289 : ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
1696 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1697 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1698 1289 : ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
1699 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1700 0 : qs_env%scf_control%outer_scf%eps_scf = qs_env%scf_control%eps_scf_hist
1701 : END IF
1702 : END IF
1703 :
1704 7757 : END SUBROUTINE wfi_set_history_variables
1705 :
1706 : ! **************************************************************************************************
1707 : !> \brief updates the snapshot buffer, taking a new snapshot
1708 : !> \param wf_history the history buffer to update
1709 : !> \param qs_env the qs_env we get the info from
1710 : !> \param dt ...
1711 : !> \par History
1712 : !> 02.2003 created [fawzi]
1713 : !> \author fawzi
1714 : ! **************************************************************************************************
1715 21299 : SUBROUTINE wfi_update(wf_history, qs_env, dt)
1716 : TYPE(qs_wf_history_type), POINTER :: wf_history
1717 : TYPE(qs_environment_type), POINTER :: qs_env
1718 : REAL(KIND=dp), INTENT(in) :: dt
1719 :
1720 21299 : CPASSERT(ASSOCIATED(wf_history))
1721 21299 : CPASSERT(wf_history%ref_count > 0)
1722 21299 : CPASSERT(ASSOCIATED(qs_env))
1723 :
1724 21299 : wf_history%snapshot_count = wf_history%snapshot_count + 1
1725 21299 : IF (wf_history%memory_depth > 0) THEN
1726 : wf_history%last_state_index = MODULO(wf_history%snapshot_count, &
1727 20512 : wf_history%memory_depth) + 1
1728 : CALL wfs_update(snapshot=wf_history%past_states &
1729 : (wf_history%last_state_index)%snapshot, wf_history=wf_history, &
1730 20512 : qs_env=qs_env, dt=dt)
1731 : END IF
1732 21299 : END SUBROUTINE wfi_update
1733 :
1734 : ! **************************************************************************************************
1735 : !> \brief reorthogonalizes the mos
1736 : !> \param qs_env the qs_env in which to orthogonalize
1737 : !> \param v_matrix the vectors to orthogonalize
1738 : !> \param n_col number of column of v to orthogonalize
1739 : !> \par History
1740 : !> 04.2003 created [fawzi]
1741 : !> \author Fawzi Mohamed
1742 : ! **************************************************************************************************
1743 28618 : SUBROUTINE reorthogonalize_vectors(qs_env, v_matrix, n_col)
1744 : TYPE(qs_environment_type), POINTER :: qs_env
1745 : TYPE(cp_fm_type), INTENT(IN) :: v_matrix
1746 : INTEGER, INTENT(in), OPTIONAL :: n_col
1747 :
1748 : CHARACTER(len=*), PARAMETER :: routineN = 'reorthogonalize_vectors'
1749 :
1750 : INTEGER :: handle, my_n_col
1751 : LOGICAL :: has_unit_metric, &
1752 : ortho_contains_cholesky, &
1753 : smearing_is_used
1754 : TYPE(cp_fm_pool_type), POINTER :: maxao_maxmo_fm_pool
1755 14309 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1756 : TYPE(dft_control_type), POINTER :: dft_control
1757 : TYPE(qs_matrix_pools_type), POINTER :: mpools
1758 : TYPE(qs_scf_env_type), POINTER :: scf_env
1759 : TYPE(scf_control_type), POINTER :: scf_control
1760 :
1761 14309 : NULLIFY (scf_env, scf_control, maxao_maxmo_fm_pool, matrix_s, mpools, dft_control)
1762 14309 : CALL timeset(routineN, handle)
1763 :
1764 14309 : CPASSERT(ASSOCIATED(qs_env))
1765 :
1766 14309 : CALL cp_fm_get_info(v_matrix, ncol_global=my_n_col)
1767 14309 : IF (PRESENT(n_col)) my_n_col = n_col
1768 : CALL get_qs_env(qs_env, mpools=mpools, &
1769 : scf_env=scf_env, &
1770 : scf_control=scf_control, &
1771 : matrix_s=matrix_s, &
1772 14309 : dft_control=dft_control)
1773 14309 : CALL mpools_get(mpools, maxao_maxmo_fm_pool=maxao_maxmo_fm_pool)
1774 14309 : IF (ASSOCIATED(scf_env)) THEN
1775 : ortho_contains_cholesky = (scf_env%method /= ot_method_nr) .AND. &
1776 : (scf_env%cholesky_method > 0) .AND. &
1777 14309 : ASSOCIATED(scf_env%ortho)
1778 : ELSE
1779 : ortho_contains_cholesky = .FALSE.
1780 : END IF
1781 :
1782 14309 : CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
1783 14309 : smearing_is_used = .FALSE.
1784 14309 : IF (dft_control%smear) THEN
1785 1146 : smearing_is_used = .TRUE.
1786 : END IF
1787 :
1788 14309 : IF (has_unit_metric) THEN
1789 3410 : CALL make_basis_simple(v_matrix, my_n_col)
1790 10899 : ELSE IF (smearing_is_used) THEN
1791 : CALL make_basis_lowdin(vmatrix=v_matrix, ncol=my_n_col, &
1792 1146 : matrix_s=matrix_s(1)%matrix)
1793 9753 : ELSE IF (ortho_contains_cholesky) THEN
1794 : CALL make_basis_cholesky(vmatrix=v_matrix, ncol=my_n_col, &
1795 6342 : ortho=scf_env%ortho)
1796 : ELSE
1797 3411 : CALL make_basis_sm(v_matrix, my_n_col, matrix_s(1)%matrix)
1798 : END IF
1799 14309 : CALL timestop(handle)
1800 14309 : END SUBROUTINE reorthogonalize_vectors
1801 :
1802 : ! **************************************************************************************************
1803 : !> \brief purges wf_history retaining only the latest snapshot
1804 : !> \param qs_env the qs env with the latest result, and that will contain
1805 : !> the purged wf_history
1806 : !> \par History
1807 : !> 05.2016 created [Nico Holmberg]
1808 : !> \author Nico Holmberg
1809 : ! **************************************************************************************************
1810 0 : SUBROUTINE wfi_purge_history(qs_env)
1811 : TYPE(qs_environment_type), POINTER :: qs_env
1812 :
1813 : CHARACTER(len=*), PARAMETER :: routineN = 'wfi_purge_history'
1814 :
1815 : INTEGER :: handle, io_unit, print_level
1816 : TYPE(cp_logger_type), POINTER :: logger
1817 : TYPE(dft_control_type), POINTER :: dft_control
1818 : TYPE(qs_wf_history_type), POINTER :: wf_history
1819 :
1820 0 : NULLIFY (dft_control, wf_history)
1821 :
1822 0 : CALL timeset(routineN, handle)
1823 0 : logger => cp_get_default_logger()
1824 0 : print_level = logger%iter_info%print_level
1825 : io_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%SCF%PRINT%PROGRAM_RUN_INFO", &
1826 0 : extension=".scfLog")
1827 :
1828 0 : CPASSERT(ASSOCIATED(qs_env))
1829 0 : CPASSERT(ASSOCIATED(qs_env%wf_history))
1830 0 : CPASSERT(qs_env%wf_history%ref_count > 0)
1831 0 : CALL get_qs_env(qs_env, dft_control=dft_control)
1832 :
1833 0 : SELECT CASE (qs_env%wf_history%interpolation_method_nr)
1834 : CASE (wfi_use_guess_method_nr, wfi_use_prev_wf_method_nr, &
1835 : wfi_use_prev_p_method_nr, wfi_use_prev_rho_r_method_nr, &
1836 : wfi_frozen_method_nr)
1837 : ! do nothing
1838 : CASE (wfi_linear_wf_method_nr, wfi_linear_p_method_nr, &
1839 : wfi_linear_ps_method_nr, wfi_ps_method_nr, &
1840 : wfi_aspc_nr, wfi_gext_proj_nr, wfi_gext_proj_qtr_nr)
1841 0 : IF (qs_env%wf_history%snapshot_count >= 2) THEN
1842 0 : IF (debug_this_module .AND. io_unit > 0) &
1843 0 : WRITE (io_unit, FMT="(T2,A)") "QS| Purging WFN history"
1844 : CALL wfi_create(wf_history, interpolation_method_nr= &
1845 : dft_control%qs_control%wf_interpolation_method_nr, &
1846 : extrapolation_order=dft_control%qs_control%wf_extrapolation_order, &
1847 0 : has_unit_metric=qs_env%has_unit_metric)
1848 : CALL set_qs_env(qs_env=qs_env, &
1849 0 : wf_history=wf_history)
1850 0 : CALL wfi_release(wf_history)
1851 0 : CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
1852 : END IF
1853 : CASE DEFAULT
1854 0 : CPABORT("Unknown extrapolation method.")
1855 : END SELECT
1856 0 : CALL timestop(handle)
1857 :
1858 0 : END SUBROUTINE wfi_purge_history
1859 :
1860 : ! **************************************************************************************************
1861 : !> \brief Gives the coefficients that best approximate the new overlap
1862 : !> as a linear combination of the previous overlaps in the
1863 : !> wf_history buffer. This is done by solving
1864 : !> argmin_a || S_{n+1} - S_{n} - \sum_i^{nvec-1} a_i (S_{n-q+i} - S_{n}) ||^2
1865 : !> \param wf_history wavefunction history buffer, containing the previous overlaps
1866 : !> \param current_overlap current overlap in dbcsr format
1867 : !> \param coeffs resulting nvec coefficients
1868 : !> \param nvec number of previous overlaps
1869 : !> \param eps Tikhonov regularization
1870 : !> \param io_unit output unit
1871 : !> \param print_level print level
1872 : !> \par History
1873 : !> 04.2026 created [Michele Nottoli]
1874 : !> \author Michele Nottoli
1875 : ! **************************************************************************************************
1876 24 : SUBROUTINE diff_fitting(wf_history, current_overlap, coeffs, nvec, eps, io_unit, print_level)
1877 : TYPE(qs_wf_history_type), POINTER :: wf_history
1878 : TYPE(dbcsr_type), INTENT(IN) :: current_overlap
1879 : INTEGER, INTENT(IN) :: nvec
1880 : REAL(KIND=dp), INTENT(OUT) :: coeffs(nvec)
1881 : REAL(KIND=dp), INTENT(IN) :: eps
1882 : INTEGER, INTENT(IN) :: io_unit, print_level
1883 :
1884 : INTEGER :: i, info, j
1885 : REAL(KIND=dp) :: error
1886 24 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: b
1887 24 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: A
1888 : TYPE(dbcsr_type) :: target_diff, tmp_i, tmp_j, tmp_k
1889 : TYPE(qs_wf_snapshot_type), POINTER :: ref_state, state
1890 :
1891 24 : IF (nvec <= 0) THEN
1892 0 : CPABORT("Not enough vectors to do the fitting")
1893 24 : ELSE IF (nvec == 1) THEN
1894 4 : coeffs(1) = 1.0_dp
1895 : RETURN
1896 : END IF
1897 :
1898 120 : ALLOCATE (A(nvec - 1, nvec - 1), b(nvec - 1))
1899 :
1900 : ! get the reference for the difference fitting
1901 20 : ref_state => wfi_get_snapshot(wf_history, wf_index=1)
1902 :
1903 : ! assemble the target difference
1904 20 : CALL dbcsr_copy(target_diff, current_overlap)
1905 20 : CALL dbcsr_add(target_diff, ref_state%overlap, 1.0_dp, -1.0_dp)
1906 :
1907 : ! allocate tmp_k
1908 20 : CALL dbcsr_copy(tmp_k, current_overlap)
1909 :
1910 : ! assemble the matrix A and the RHS b
1911 76 : DO i = 2, nvec
1912 56 : state => wfi_get_snapshot(wf_history, wf_index=i)
1913 56 : CALL dbcsr_copy(tmp_i, state%overlap)
1914 56 : CALL dbcsr_add(tmp_i, ref_state%overlap, 1.0_dp, -1.0_dp)
1915 56 : CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_i, target_diff, 0.0_dp, tmp_k)
1916 56 : CALL dbcsr_trace(tmp_k, b(i - 1))
1917 :
1918 196 : DO j = 2, i
1919 120 : state => wfi_get_snapshot(wf_history, wf_index=j)
1920 120 : CALL dbcsr_copy(tmp_j, state%overlap)
1921 120 : CALL dbcsr_add(tmp_j, ref_state%overlap, 1.0_dp, -1.0_dp)
1922 120 : CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_i, tmp_j, 0.0_dp, tmp_k)
1923 120 : CALL dbcsr_trace(tmp_k, A(j - 1, i - 1))
1924 176 : A(i - 1, j - 1) = A(j - 1, i - 1)
1925 : END DO
1926 : END DO
1927 :
1928 : ! add the Tikhonov regularization
1929 76 : DO i = 1, nvec - 1
1930 76 : A(i, i) = A(i, i) + eps**2
1931 : END DO
1932 :
1933 : ! solve the linear system
1934 20 : CALL dposv('u', nvec - 1, 1, A, nvec - 1, b, nvec - 1, info)
1935 20 : IF (info /= 0) THEN
1936 0 : CPABORT("DPOSV failed.")
1937 : END IF
1938 :
1939 : ! set the coefficient for the reference snapshot
1940 76 : coeffs(1) = 1.0_dp - SUM(b)
1941 76 : coeffs(2:nvec) = b(:)
1942 :
1943 : ! as a consistency check, print how well the current overlap
1944 : ! is approximated by the linear combination of previous overlaps
1945 20 : IF (print_level > low_print_level) THEN
1946 20 : CALL dbcsr_copy(tmp_i, current_overlap)
1947 96 : DO i = 1, nvec
1948 76 : state => wfi_get_snapshot(wf_history, wf_index=i)
1949 96 : CALL dbcsr_add(tmp_i, state%overlap, 1.0_dp, -coeffs(i))
1950 : END DO
1951 20 : error = dbcsr_frobenius_norm(tmp_i)/dbcsr_frobenius_norm(state%overlap)
1952 20 : IF (io_unit > 0) THEN
1953 10 : WRITE (UNIT=io_unit, FMT="(/,T2,A,F20.10)") "GEXT overlap fitting error:", error
1954 : END IF
1955 : END IF
1956 :
1957 : ! free the memory
1958 20 : CALL dbcsr_release(tmp_i)
1959 20 : CALL dbcsr_release(tmp_j)
1960 20 : CALL dbcsr_release(tmp_k)
1961 20 : CALL dbcsr_release(target_diff)
1962 20 : DEALLOCATE (A, b)
1963 :
1964 24 : END SUBROUTINE diff_fitting
1965 :
1966 : ! **************************************************************************************************
1967 : !> \brief Gives the coefficients that best approximate the new overlap
1968 : !> as a time reversible linear combination of the previous overlaps in the
1969 : !> wf_history buffer. This is done by solving
1970 : !> argmin_a || S_{n+1} + S_{n+1-nvec}
1971 : !> - \sum_{i=1}^q a_i (S_{n+1-nvec+i} + S_{n+1-i}) ||^2
1972 : !> with q = nvec/2 if nvec is even, or q = (nvec-1)/2 if odd.
1973 : !> \param wf_history wavefunction history buffer, containing the previous overlaps
1974 : !> \param current_overlap current overlap in dbcsr format
1975 : !> \param coeffs resulting nvec coefficients
1976 : !> \param nvec number of previous overlaps
1977 : !> \param eps Tikhonov regularization
1978 : !> \param io_unit output unit
1979 : !> \param print_level print level
1980 : !> \par History
1981 : !> 04.2026 created [Michele Nottoli]
1982 : ! **************************************************************************************************
1983 24 : SUBROUTINE tr_fitting(wf_history, current_overlap, coeffs, nvec, eps, io_unit, print_level)
1984 : TYPE(qs_wf_history_type), POINTER :: wf_history
1985 : TYPE(dbcsr_type), INTENT(IN) :: current_overlap
1986 : INTEGER, INTENT(IN) :: nvec
1987 : REAL(KIND=dp), INTENT(OUT) :: coeffs(nvec)
1988 : REAL(KIND=dp), INTENT(IN) :: eps
1989 : INTEGER, INTENT(IN) :: io_unit, print_level
1990 :
1991 : INTEGER :: i, info, j, ntr
1992 : REAL(KIND=dp) :: error
1993 24 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: b
1994 24 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: A
1995 : TYPE(dbcsr_type) :: target_overlap, tmp_i, tmp_j, tmp_k
1996 : TYPE(qs_wf_snapshot_type), POINTER :: ref_state, state
1997 :
1998 24 : IF (nvec <= 0) THEN
1999 0 : CPABORT("Not enough vectors to do the fitting")
2000 24 : ELSE IF (nvec == 1) THEN
2001 4 : coeffs(1) = 1.0_dp
2002 : RETURN
2003 : END IF
2004 :
2005 20 : IF (MOD(nvec, 2) == 0) THEN
2006 8 : ntr = nvec/2
2007 : ELSE
2008 12 : ntr = (nvec - 1)/2
2009 : END IF
2010 :
2011 120 : ALLOCATE (A(ntr, ntr), b(ntr))
2012 :
2013 : ! get the reference for the difference fitting
2014 20 : ref_state => wfi_get_snapshot(wf_history, wf_index=nvec)
2015 :
2016 : ! assemble the target sum
2017 20 : CALL dbcsr_copy(target_overlap, current_overlap)
2018 20 : CALL dbcsr_add(target_overlap, ref_state%overlap, 1.0_dp, 1.0_dp)
2019 :
2020 : ! allocate tmp_k
2021 20 : CALL dbcsr_copy(tmp_k, current_overlap)
2022 :
2023 : ! assemble the matrix A and the RHS b
2024 52 : DO i = 1, ntr
2025 32 : state => wfi_get_snapshot(wf_history, wf_index=i)
2026 32 : CALL dbcsr_copy(tmp_i, state%overlap)
2027 32 : state => wfi_get_snapshot(wf_history, wf_index=nvec - i)
2028 32 : CALL dbcsr_add(tmp_i, state%overlap, 1.0_dp, 1.0_dp)
2029 :
2030 : CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_i, target_overlap, &
2031 32 : 0.0_dp, tmp_k)
2032 32 : CALL dbcsr_trace(tmp_k, b(i))
2033 96 : DO j = 1, i
2034 44 : state => wfi_get_snapshot(wf_history, wf_index=j)
2035 44 : CALL dbcsr_copy(tmp_j, state%overlap)
2036 44 : state => wfi_get_snapshot(wf_history, wf_index=nvec - j)
2037 44 : CALL dbcsr_add(tmp_j, state%overlap, 1.0_dp, 1.0_dp)
2038 44 : CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_i, tmp_j, 0.0_dp, tmp_k)
2039 44 : CALL dbcsr_trace(tmp_k, A(j, i))
2040 76 : A(i, j) = A(j, i)
2041 : END DO
2042 : END DO
2043 :
2044 : ! add the Tikhonov regularization
2045 52 : DO i = 1, ntr
2046 52 : A(i, i) = A(i, i) + eps**2
2047 : END DO
2048 :
2049 : ! solve the linear system
2050 20 : CALL dposv('u', ntr, 1, A, ntr, b, ntr, info)
2051 20 : IF (info /= 0) THEN
2052 0 : CPABORT("DPOSV failed.")
2053 : END IF
2054 :
2055 : ! reorder the coefficients
2056 96 : coeffs = 0.0_dp
2057 20 : coeffs(nvec) = -1.0_dp
2058 52 : DO i = 1, ntr
2059 32 : coeffs(i) = coeffs(i) + b(i)
2060 52 : coeffs(nvec - i) = coeffs(nvec - i) + b(i)
2061 : END DO
2062 :
2063 : ! as a consistency check, print how well the current overlap
2064 : ! is approximated by the linear combination of previous overlaps
2065 20 : IF (print_level > low_print_level) THEN
2066 20 : CALL dbcsr_copy(tmp_i, current_overlap)
2067 96 : DO i = 1, nvec
2068 76 : state => wfi_get_snapshot(wf_history, wf_index=i)
2069 96 : CALL dbcsr_add(tmp_i, state%overlap, 1.0_dp, -coeffs(i))
2070 : END DO
2071 20 : error = dbcsr_frobenius_norm(tmp_i)/dbcsr_frobenius_norm(state%overlap)
2072 20 : IF (io_unit > 0) THEN
2073 10 : WRITE (UNIT=io_unit, FMT="(/,T2,A,F20.10)") "GEXT overlap fitting error:", error
2074 : END IF
2075 : END IF
2076 :
2077 : ! free the memory
2078 20 : CALL dbcsr_release(tmp_i)
2079 20 : CALL dbcsr_release(tmp_j)
2080 20 : CALL dbcsr_release(tmp_k)
2081 20 : CALL dbcsr_release(target_overlap)
2082 20 : DEALLOCATE (A, b)
2083 :
2084 24 : END SUBROUTINE tr_fitting
2085 :
2086 : END MODULE qs_wf_history_methods
|
