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