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