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