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