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