(git:85b8a9b)
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)) &
765 extrapolation_method_nr = actual_extrapolation_method_nr
766 my_orthogonal_wf = .false.
767
768 SELECT CASE (actual_extrapolation_method_nr)
769 CASE (wfi_frozen_method_nr)
770 cpassert(.NOT. do_kpoints)
771 t0_state => wfi_get_snapshot(wf_history, wf_index=1)
772 cpassert(ASSOCIATED(t0_state%rho_frozen))
773
774 nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
775 CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
776
777 CALL qs_rho_get(t0_state%rho_frozen, rho_ao=rho_frozen_ao)
778 CALL qs_rho_get(rho, rho_ao=rho_ao)
779 DO ispin = 1, SIZE(rho_frozen_ao)
780 CALL dbcsr_copy(rho_ao(ispin)%matrix, &
781 rho_frozen_ao(ispin)%matrix, &
782 keep_sparsity=.true.)
783 END DO
784 !FM updating rho_ao directly with t0_state%rho_ao would have the
785 !FM wrong matrix structure
786 CALL qs_rho_update_rho(rho, qs_env=qs_env)
787 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
788
789 my_orthogonal_wf = .false.
790 CASE (wfi_use_prev_rho_r_method_nr, wfi_use_prev_p_method_nr)
791 IF (actual_extrapolation_method_nr == wfi_use_prev_rho_r_method_nr) THEN
792 CALL cp_warn(__location__, &
793 "USE_PREV_RHO_R is deprecated and will be removed in a future "// &
794 "release; it is now an alias of USE_PREV_P.")
795 END IF
796 t0_state => wfi_get_snapshot(wf_history, wf_index=1)
797 nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
798 CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
799 IF (do_kpoints) THEN
800 cpassert(ASSOCIATED(t0_state%rho_ao_kp))
801 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
802 DO ispin = 1, SIZE(t0_state%rho_ao_kp, 1)
803 DO img = 1, SIZE(t0_state%rho_ao_kp, 2)
804 IF (img > SIZE(rho_ao_kp, 2)) THEN
805 cpwarn("Change in cell neighborlist: might affect quality of initial guess")
806 ELSE
807 CALL dbcsr_copy(rho_ao_kp(ispin, img)%matrix, &
808 t0_state%rho_ao_kp(ispin, img)%matrix, &
809 keep_sparsity=.true.)
810 END IF
811 END DO
812 END DO
813 ELSE
814 cpassert(ASSOCIATED(t0_state%rho_ao))
815 CALL qs_rho_get(rho, rho_ao=rho_ao)
816 DO ispin = 1, SIZE(t0_state%rho_ao)
817 CALL dbcsr_copy(rho_ao(ispin)%matrix, &
818 t0_state%rho_ao(ispin)%matrix, &
819 keep_sparsity=.true.)
820 END DO
821 END IF
822 !FM updating rho_ao directly with t0_state%rho_ao would have the
823 !FM wrong matrix structure
824 CALL qs_rho_update_rho(rho, qs_env=qs_env)
825 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
826 CASE (wfi_use_prev_wf_method_nr)
827 my_orthogonal_wf = .true.
828 nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
829 CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
830
831 IF (do_kpoints) THEN
832 CALL wfi_use_prev_wf_kp(qs_env, io_unit, print_level)
833 ELSE
834 CALL qs_rho_get(rho, rho_ao=rho_ao)
835 DO ispin = 1, SIZE(mos)
836 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
837 CALL reorthogonalize_vectors(qs_env, v_matrix=mo_coeff, n_col=nmo)
838 CALL calculate_density_matrix(mo_set=mos(ispin), density_matrix=rho_ao(ispin)%matrix)
839 END DO
840 CALL qs_rho_update_rho(rho, qs_env=qs_env)
841 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
842 END IF
843
844 CASE (wfi_use_guess_method_nr)
845 !FM more clean to do it here, but it
846 !FM might need to read a file (restart) and thus globenv
847 !FM I do not want globenv here, thus done by the caller
848 !FM (btw. it also needs the eigensolver, and unless you relocate it
849 !FM gives circular dependencies)
850 nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
851 CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
852 CASE (wfi_linear_wf_method_nr)
853 cpassert(.NOT. do_kpoints)
854 t0_state => wfi_get_snapshot(wf_history, wf_index=2)
855 t1_state => wfi_get_snapshot(wf_history, wf_index=1)
856 cpassert(ASSOCIATED(t0_state))
857 cpassert(ASSOCIATED(t1_state))
858 cpassert(ASSOCIATED(t0_state%wf))
859 cpassert(ASSOCIATED(t1_state%wf))
860 nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
861 CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
862
863 my_orthogonal_wf = .true.
864 t0 = 0.0_dp
865 t1 = t1_state%dt
866 t2 = t1 + dt
867 CALL qs_rho_get(rho, rho_ao=rho_ao)
868 DO ispin = 1, SIZE(mos)
869 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, &
870 nmo=nmo)
871 CALL cp_fm_scale_and_add(alpha=0.0_dp, &
872 matrix_a=mo_coeff, &
873 matrix_b=t1_state%wf(ispin), &
874 beta=(t2 - t0)/(t1 - t0))
875 ! this copy should be unnecessary
876 CALL cp_fm_scale_and_add(alpha=1.0_dp, &
877 matrix_a=mo_coeff, &
878 beta=(t1 - t2)/(t1 - t0), matrix_b=t0_state%wf(ispin))
879 CALL reorthogonalize_vectors(qs_env, &
880 v_matrix=mo_coeff, &
881 n_col=nmo)
882 CALL calculate_density_matrix(mo_set=mos(ispin), &
883 density_matrix=rho_ao(ispin)%matrix)
884 END DO
885 CALL qs_rho_update_rho(rho, qs_env=qs_env)
886
887 CALL qs_ks_did_change(qs_env%ks_env, &
888 rho_changed=.true.)
889 CASE (wfi_linear_p_method_nr)
890 t0_state => wfi_get_snapshot(wf_history, wf_index=2)
891 t1_state => wfi_get_snapshot(wf_history, wf_index=1)
892 cpassert(ASSOCIATED(t0_state))
893 cpassert(ASSOCIATED(t1_state))
894 IF (do_kpoints) THEN
895 cpassert(ASSOCIATED(t0_state%rho_ao_kp))
896 cpassert(ASSOCIATED(t1_state%rho_ao_kp))
897 ELSE
898 cpassert(ASSOCIATED(t0_state%rho_ao))
899 cpassert(ASSOCIATED(t1_state%rho_ao))
900 END IF
901 nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
902 CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
903
904 t0 = 0.0_dp
905 t1 = t1_state%dt
906 t2 = t1 + dt
907 IF (do_kpoints) THEN
908 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
909 DO ispin = 1, SIZE(rho_ao_kp, 1)
910 DO img = 1, SIZE(rho_ao_kp, 2)
911 IF (img > SIZE(t0_state%rho_ao_kp, 2) .OR. &
912 img > SIZE(t1_state%rho_ao_kp, 2)) THEN
913 cpwarn("Change in cell neighborlist: might affect quality of initial guess")
914 ELSE
915 CALL dbcsr_add(rho_ao_kp(ispin, img)%matrix, t1_state%rho_ao_kp(ispin, img)%matrix, &
916 alpha_scalar=0.0_dp, beta_scalar=(t2 - t0)/(t1 - t0)) ! this copy should be unnecessary
917 CALL dbcsr_add(rho_ao_kp(ispin, img)%matrix, t0_state%rho_ao_kp(ispin, img)%matrix, &
918 alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0))
919 END IF
920 END DO
921 END DO
922 ELSE
923 CALL qs_rho_get(rho, rho_ao=rho_ao)
924 DO ispin = 1, SIZE(rho_ao)
925 CALL dbcsr_add(rho_ao(ispin)%matrix, t1_state%rho_ao(ispin)%matrix, &
926 alpha_scalar=0.0_dp, beta_scalar=(t2 - t0)/(t1 - t0)) ! this copy should be unnecessary
927 CALL dbcsr_add(rho_ao(ispin)%matrix, t0_state%rho_ao(ispin)%matrix, &
928 alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0))
929 END DO
930 END IF
931 CALL qs_rho_update_rho(rho, qs_env=qs_env)
932 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
933 CASE (wfi_linear_ps_method_nr)
934 ! wf not calculated, extract with PSC renormalized?
935 ! use wf_linear?
936 cpassert(.NOT. do_kpoints)
937 t0_state => wfi_get_snapshot(wf_history, wf_index=2)
938 t1_state => wfi_get_snapshot(wf_history, wf_index=1)
939 cpassert(ASSOCIATED(t0_state))
940 cpassert(ASSOCIATED(t1_state))
941 cpassert(ASSOCIATED(t0_state%wf))
942 cpassert(ASSOCIATED(t1_state%wf))
943 IF (wf_history%store_overlap) THEN
944 cpassert(ASSOCIATED(t0_state%overlap))
945 cpassert(ASSOCIATED(t1_state%overlap))
946 END IF
947 nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
948 IF (nvec >= wf_history%memory_depth) THEN
949 IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. (qs_env%scf_control%eps_scf_hist /= 0)) THEN
950 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
951 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
952 qs_env%scf_control%outer_scf%have_scf = .false.
953 ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
954 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
955 qs_env%scf_control%outer_scf%have_scf = .false.
956 ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
957 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
958 END IF
959 END IF
960
961 my_orthogonal_wf = .true.
962 ! use PS_2=2 PS_1-PS_0
963 ! C_2 comes from using PS_2 as a projector acting on C_1
964 CALL qs_rho_get(rho, rho_ao=rho_ao)
965 DO ispin = 1, SIZE(mos)
966 NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
967 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
968 CALL cp_fm_get_info(mo_coeff, nrow_global=n, ncol_global=k, &
969 matrix_struct=matrix_struct)
970 CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, &
971 nrow_global=k, ncol_global=k)
972 CALL cp_fm_create(csc, matrix_struct_new)
973 CALL cp_fm_struct_release(matrix_struct_new)
974
975 IF (use_overlap) THEN
976 CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), mo_coeff, k)
977 CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), mo_coeff, 0.0_dp, csc)
978 ELSE
979 CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
980 t1_state%wf(ispin), 0.0_dp, csc)
981 END IF
982 CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, mo_coeff)
983 CALL cp_fm_release(csc)
984 CALL cp_fm_scale_and_add(-1.0_dp, mo_coeff, 2.0_dp, t1_state%wf(ispin))
985 CALL reorthogonalize_vectors(qs_env, &
986 v_matrix=mo_coeff, &
987 n_col=k)
988 CALL calculate_density_matrix(mo_set=mos(ispin), &
989 density_matrix=rho_ao(ispin)%matrix)
990 END DO
991 CALL qs_rho_update_rho(rho, qs_env=qs_env)
992 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
993
994 CASE (wfi_ps_method_nr)
995 ! figure out the actual number of vectors to use in the extrapolation:
996 nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
997 cpassert(nvec > 0)
998 IF (nvec >= wf_history%memory_depth) THEN
999 IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. (qs_env%scf_control%eps_scf_hist /= 0)) THEN
1000 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1001 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1002 qs_env%scf_control%outer_scf%have_scf = .false.
1003 ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
1004 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1005 qs_env%scf_control%outer_scf%have_scf = .false.
1006 ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
1007 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1008 END IF
1009 END IF
1010
1011 IF (do_kpoints) THEN
1012 CALL wfi_extrapolate_ps_aspc_kp(wf_history, qs_env, nvec, io_unit, print_level)
1013 my_orthogonal_wf = .true.
1014 ELSE
1015 my_orthogonal_wf = .true.
1016 DO ispin = 1, SIZE(mos)
1017 NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
1018 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1019 CALL cp_fm_get_info(mo_coeff, nrow_global=n, ncol_global=k, &
1020 matrix_struct=matrix_struct)
1021 CALL cp_fm_create(fm_tmp, matrix_struct)
1022 CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, &
1023 nrow_global=k, ncol_global=k)
1024 CALL cp_fm_create(csc, matrix_struct_new)
1025 CALL cp_fm_struct_release(matrix_struct_new)
1026 ! first the most recent
1027 t1_state => wfi_get_snapshot(wf_history, wf_index=1)
1028 CALL cp_fm_to_fm(t1_state%wf(ispin), mo_coeff)
1029 alpha = nvec
1030 CALL cp_fm_scale(alpha, mo_coeff)
1031 CALL qs_rho_get(rho, rho_ao=rho_ao)
1032 DO i = 2, nvec
1033 t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1034 IF (use_overlap) THEN
1035 CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
1036 CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
1037 ELSE
1038 CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
1039 t1_state%wf(ispin), 0.0_dp, csc)
1040 END IF
1041 CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
1042 alpha = -1.0_dp*alpha*real(nvec - i + 1, dp)/real(i, dp)
1043 CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp)
1044 END DO
1045
1046 CALL cp_fm_release(csc)
1047 CALL cp_fm_release(fm_tmp)
1048 CALL reorthogonalize_vectors(qs_env, &
1049 v_matrix=mo_coeff, &
1050 n_col=k)
1051 CALL calculate_density_matrix(mo_set=mos(ispin), &
1052 density_matrix=rho_ao(ispin)%matrix)
1053 END DO
1054 CALL qs_rho_update_rho(rho, qs_env=qs_env)
1055 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
1056 END IF
1057
1058 CASE (wfi_aspc_nr)
1059 CALL cite_reference(kolafa2004)
1060 CALL cite_reference(kuhne2007)
1061 ! figure out the actual number of vectors to use in the extrapolation:
1062 nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
1063 cpassert(nvec > 0)
1064 IF (nvec >= wf_history%memory_depth) THEN
1065 IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. &
1066 (qs_env%scf_control%eps_scf_hist /= 0)) THEN
1067 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1068 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1069 qs_env%scf_control%outer_scf%have_scf = .false.
1070 ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
1071 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1072 qs_env%scf_control%outer_scf%have_scf = .false.
1073 ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
1074 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1075 END IF
1076 END IF
1077
1078 IF (do_kpoints) THEN
1079 CALL wfi_extrapolate_ps_aspc_kp(wf_history, qs_env, nvec, io_unit, print_level)
1080 my_orthogonal_wf = .true.
1081 ELSE
1082 my_orthogonal_wf = .true.
1083 CALL qs_rho_get(rho, rho_ao=rho_ao)
1084 DO ispin = 1, SIZE(mos)
1085 NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
1086 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1087 CALL cp_fm_get_info(mo_coeff, &
1088 nrow_global=n, &
1089 ncol_global=k, &
1090 matrix_struct=matrix_struct)
1091 CALL cp_fm_create(fm_tmp, matrix_struct, set_zero=.true.)
1092 CALL cp_fm_struct_create(matrix_struct_new, &
1093 template_fmstruct=matrix_struct, &
1094 nrow_global=k, &
1095 ncol_global=k)
1096 CALL cp_fm_create(csc, matrix_struct_new, set_zero=.true.)
1097 CALL cp_fm_struct_release(matrix_struct_new)
1098 ! first the most recent
1099 t1_state => wfi_get_snapshot(wf_history, &
1100 wf_index=1)
1101 CALL cp_fm_to_fm(t1_state%wf(ispin), mo_coeff)
1102 alpha = real(4*nvec - 2, kind=dp)/real(nvec + 1, kind=dp)
1103 IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1104 WRITE (unit=io_unit, fmt="(/,T2,A,/,/,T3,A,I0,/,/,T3,A2,I0,A4,F10.6)") &
1105 "Parameters for the always stable predictor-corrector (ASPC) method:", &
1106 "ASPC order: ", max(nvec - 2, 0), &
1107 "B(", 1, ") = ", alpha
1108 END IF
1109 CALL cp_fm_scale(alpha, mo_coeff)
1110
1111 DO i = 2, nvec
1112 t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1113 IF (use_overlap) THEN
1114 CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
1115 CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
1116 ELSE
1117 CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
1118 t1_state%wf(ispin), 0.0_dp, csc)
1119 END IF
1120 CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
1121 alpha = (-1.0_dp)**(i + 1)*real(i, kind=dp)* &
1122 binomial(2*nvec, nvec - i)/binomial(2*nvec - 2, nvec - 1)
1123 IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1124 WRITE (unit=io_unit, fmt="(T3,A2,I0,A4,F10.6)") &
1125 "B(", i, ") = ", alpha
1126 END IF
1127 CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp)
1128 END DO
1129 CALL cp_fm_release(csc)
1130 CALL cp_fm_release(fm_tmp)
1131 CALL reorthogonalize_vectors(qs_env, &
1132 v_matrix=mo_coeff, &
1133 n_col=k)
1134 CALL calculate_density_matrix(mo_set=mos(ispin), &
1135 density_matrix=rho_ao(ispin)%matrix)
1136 END DO
1137 CALL qs_rho_update_rho(rho, qs_env=qs_env)
1138 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
1139 END IF ! do_kpoints
1140
1141 CASE (wfi_gext_proj_nr)
1142 IF (do_kpoints) THEN
1143 nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
1144 cpassert(nvec > 0)
1145 CALL wfi_extrapolate_gext_proj_kp(wf_history, qs_env, nvec, io_unit, print_level)
1146 ELSE
1147
1148 ! figure out the actual number of vectors to use in the extrapolation:
1149 nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
1150 IF (nvec >= wf_history%memory_depth) THEN
1151 IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. &
1152 (qs_env%scf_control%eps_scf_hist /= 0)) THEN
1153 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1154 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1155 qs_env%scf_control%outer_scf%have_scf = .false.
1156 ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
1157 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1158 qs_env%scf_control%outer_scf%have_scf = .false.
1159 ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
1160 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1161 END IF
1162 END IF
1163 cpassert(nvec > 0)
1164
1165 ! get the coefficients for the fitting
1166 ALLOCATE (coeffs(nvec))
1167 NULLIFY (matrix_s)
1168 CALL get_qs_env(qs_env, matrix_s=matrix_s)
1169 CALL diff_fitting(wf_history, matrix_s(1)%matrix, coeffs, nvec, &
1170 1e-4_dp, io_unit, print_level)
1171
1172 my_orthogonal_wf = .true.
1173 CALL qs_rho_get(rho, rho_ao=rho_ao)
1174 DO ispin = 1, SIZE(mos)
1175 NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
1176 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1177 CALL cp_fm_get_info(mo_coeff, &
1178 nrow_global=n, &
1179 ncol_global=k, &
1180 matrix_struct=matrix_struct)
1181 CALL cp_fm_create(fm_tmp, matrix_struct)
1182 CALL cp_fm_struct_create(matrix_struct_new, &
1183 template_fmstruct=matrix_struct, &
1184 nrow_global=k, &
1185 ncol_global=k)
1186 CALL cp_fm_create(csc, matrix_struct_new)
1187 CALL cp_fm_struct_release(matrix_struct_new)
1188
1189 t1_state => wfi_get_snapshot(wf_history, wf_index=1)
1190
1191 ! do the linear combination of previous PSs
1192 CALL cp_fm_set_all(mo_coeff, 0.0_dp)
1193 DO i = 1, nvec
1194 t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1195 CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
1196 CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
1197 CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
1198 CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, coeffs(i), fm_tmp)
1199 END DO
1200 CALL cp_fm_release(csc)
1201 CALL cp_fm_release(fm_tmp)
1202 CALL reorthogonalize_vectors(qs_env, &
1203 v_matrix=mo_coeff, &
1204 n_col=k)
1205 CALL calculate_density_matrix(mo_set=mos(ispin), &
1206 density_matrix=rho_ao(ispin)%matrix)
1207 END DO
1208 CALL qs_rho_update_rho(rho, qs_env=qs_env)
1209 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
1210
1211 DEALLOCATE (coeffs)
1212
1213 END IF
1214
1215 CASE (wfi_gext_proj_qtr_nr)
1216 IF (do_kpoints) THEN
1217 nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
1218 cpassert(nvec > 0)
1219 CALL wfi_extrapolate_gext_proj_kp(wf_history, qs_env, nvec, io_unit, print_level)
1220 ELSE
1221
1222 ! figure out the actual number of vectors to use in the extrapolation:
1223 nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
1224 IF (nvec >= wf_history%memory_depth) THEN
1225 IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. &
1226 (qs_env%scf_control%eps_scf_hist /= 0)) THEN
1227 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1228 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1229 qs_env%scf_control%outer_scf%have_scf = .false.
1230 ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
1231 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1232 qs_env%scf_control%outer_scf%have_scf = .false.
1233 ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
1234 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1235 END IF
1236 END IF
1237 cpassert(nvec > 0)
1238
1239 ! get the coefficients for the fitting
1240 ALLOCATE (coeffs(nvec))
1241 NULLIFY (matrix_s)
1242 CALL get_qs_env(qs_env, matrix_s=matrix_s)
1243 CALL tr_fitting(wf_history, matrix_s(1)%matrix, coeffs, nvec, &
1244 1e-4_dp, io_unit, print_level)
1245
1246 my_orthogonal_wf = .true.
1247 CALL qs_rho_get(rho, rho_ao=rho_ao)
1248 DO ispin = 1, SIZE(mos)
1249 NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
1250 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1251 CALL cp_fm_get_info(mo_coeff, &
1252 nrow_global=n, &
1253 ncol_global=k, &
1254 matrix_struct=matrix_struct)
1255 CALL cp_fm_create(fm_tmp, matrix_struct)
1256 CALL cp_fm_struct_create(matrix_struct_new, &
1257 template_fmstruct=matrix_struct, &
1258 nrow_global=k, &
1259 ncol_global=k)
1260 CALL cp_fm_create(csc, matrix_struct_new)
1261 CALL cp_fm_struct_release(matrix_struct_new)
1262
1263 t1_state => wfi_get_snapshot(wf_history, wf_index=1)
1264
1265 ! do the linear combination of previous PSs
1266 CALL cp_fm_set_all(mo_coeff, 0.0_dp)
1267 DO i = 1, nvec
1268 t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1269 CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
1270 CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
1271 CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
1272 CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, coeffs(i), fm_tmp)
1273 END DO
1274 CALL cp_fm_release(csc)
1275 CALL cp_fm_release(fm_tmp)
1276 CALL reorthogonalize_vectors(qs_env, &
1277 v_matrix=mo_coeff, &
1278 n_col=k)
1279 CALL calculate_density_matrix(mo_set=mos(ispin), &
1280 density_matrix=rho_ao(ispin)%matrix)
1281 END DO
1282 CALL qs_rho_update_rho(rho, qs_env=qs_env)
1283 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
1284
1285 DEALLOCATE (coeffs)
1286
1287 END IF
1288
1289 CASE default
1290 CALL cp_abort(__location__, &
1291 "Unknown interpolation method: "// &
1292 trim(adjustl(cp_to_string(wf_history%interpolation_method_nr))))
1293 END SELECT
1294 IF (PRESENT(orthogonal_wf)) orthogonal_wf = my_orthogonal_wf
1295 CALL cp_print_key_finished_output(io_unit, logger, qs_env%input, &
1296 "DFT%SCF%PRINT%PROGRAM_RUN_INFO")
1297 CALL timestop(handle)
1298 END SUBROUTINE wfi_extrapolate
1299
1300! **************************************************************************************************
1301!> \brief Reorthogonalizes the wavefunctions from the previous step for k-points
1302!> using the current S(k) metric and rebuilds the density matrix.
1303!> \param qs_env The QS environment
1304!> \param io_unit output unit
1305!> \param print_level print level
1306!> \param pbc_shift_ref ...
1307!> \param load_snapshot_wf ...
1308! **************************************************************************************************
1309 SUBROUTINE wfi_use_prev_wf_kp(qs_env, io_unit, print_level, pbc_shift_ref, load_snapshot_wf)
1310 TYPE(qs_environment_type), POINTER :: qs_env
1311 INTEGER, INTENT(IN) :: io_unit, print_level
1312 INTEGER, DIMENSION(:, :), INTENT(IN), OPTIONAL :: pbc_shift_ref
1313 LOGICAL, INTENT(IN), OPTIONAL :: load_snapshot_wf
1314
1315 CHARACTER(len=*), PARAMETER :: routinen = 'wfi_use_prev_wf_kp'
1316
1317 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: col_scaling
1318 INTEGER :: chol_info, handle, igroup, ik, ikp, &
1319 indx, ispin, j, kplocal, nao, nkp, &
1320 nkp_groups, nmo, nspin
1321 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: pbc_shift_cur, pbc_shift_src
1322 INTEGER, DIMENSION(2) :: kp_range
1323 INTEGER, DIMENSION(:, :), POINTER :: kp_dist
1324 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1325 LOGICAL :: my_kpgrp, reload_snapshot_wf, &
1326 use_pbc_phase_ref, use_real_wfn
1327 REAL(kind=dp) :: eval_thresh
1328 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1329 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
1330 TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
1331 TYPE(cp_cfm_type) :: cfm_evecs, cfm_mhalf, cfm_nao_nmo_work, &
1332 cmos_new, csc_cfm
1333 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: csmat_cur
1334 TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools_kp
1335 TYPE(cp_fm_struct_type), POINTER :: ao_ao_struct, nmo_nmo_struct
1336 TYPE(cp_fm_type) :: fmdummy, fmlocal
1337 TYPE(cp_fm_type), POINTER :: imos, mo_coeff, rmos
1338 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, rho_ao_kp
1339 TYPE(dbcsr_type), POINTER :: cmatrix_db, rmatrix, tmpmat
1340 TYPE(dft_control_type), POINTER :: dft_control
1341 TYPE(kpoint_env_type), POINTER :: kp
1342 TYPE(kpoint_type), POINTER :: kpoints
1343 TYPE(mp_para_env_type), POINTER :: para_env
1344 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1345 POINTER :: sab_nl
1346 TYPE(qs_matrix_pools_type), POINTER :: mpools_kp
1347 TYPE(qs_rho_type), POINTER :: rho
1348 TYPE(qs_scf_env_type), POINTER :: scf_env
1349 TYPE(qs_wf_history_type), POINTER :: wf_history
1350 TYPE(qs_wf_snapshot_type), POINTER :: t1_state
1351 TYPE(scf_control_type), POINTER :: scf_control
1352
1353 CALL timeset(routinen, handle)
1354
1355 NULLIFY (kpoints, dft_control, matrix_s_kp, scf_env, scf_control, rho, sab_nl, kp, &
1356 mo_coeff, rmos, imos, wf_history, t1_state)
1357
1358 CALL get_qs_env(qs_env, kpoints=kpoints, dft_control=dft_control, &
1359 matrix_s_kp=matrix_s_kp, scf_env=scf_env, &
1360 scf_control=scf_control, rho=rho)
1361 CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn, &
1362 kp_range=kp_range, nkp_groups=nkp_groups, kp_dist=kp_dist, &
1363 sab_nl=sab_nl, cell_to_index=cell_to_index)
1364 kplocal = kp_range(2) - kp_range(1) + 1
1365
1366 IF (use_real_wfn) THEN
1367 CALL timestop(handle)
1368 RETURN
1369 END IF
1370
1371 wf_history => qs_env%wf_history
1372 reload_snapshot_wf = .false.
1373 IF (PRESENT(load_snapshot_wf)) reload_snapshot_wf = load_snapshot_wf
1374 IF (PRESENT(pbc_shift_ref)) THEN
1375 ALLOCATE (pbc_shift_src(3, SIZE(pbc_shift_ref, 2)))
1376 pbc_shift_src(:, :) = pbc_shift_ref(:, :)
1377 use_pbc_phase_ref = .true.
1378 ELSE
1379 use_pbc_phase_ref = .false.
1380 IF (ASSOCIATED(wf_history)) THEN
1381 IF (wf_history%store_wf_kp .AND. wf_history%snapshot_count > 0) THEN
1382 t1_state => wfi_get_snapshot(wf_history, wf_index=1)
1383 cpassert(ASSOCIATED(t1_state%wf_kp))
1384 cpassert(ASSOCIATED(t1_state%kp_pbc_shift))
1385 reload_snapshot_wf = .true.
1386 ALLOCATE (pbc_shift_src(3, SIZE(t1_state%kp_pbc_shift, 2)))
1387 pbc_shift_src(:, :) = t1_state%kp_pbc_shift(:, :)
1388 use_pbc_phase_ref = .true.
1389 END IF
1390 END IF
1391 END IF
1392 IF (use_pbc_phase_ref) CALL wfi_compute_kp_pbc_shift(qs_env, pbc_shift_cur)
1393
1394 kp => kpoints%kp_env(1)%kpoint_env
1395 nspin = SIZE(kp%mos, 2)
1396 CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo, mo_coeff=mo_coeff)
1397
1398 IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1399 WRITE (unit=io_unit, fmt="(/,T2,A)") &
1400 "Using previous wavefunctions as initial guess for k-points (with reorthogonalization)"
1401 END IF
1402
1403 ! Allocate dbcsr work matrices
1404 ALLOCATE (rmatrix, cmatrix_db, tmpmat)
1405 CALL dbcsr_create(rmatrix, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
1406 CALL dbcsr_create(cmatrix_db, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
1407 CALL dbcsr_create(tmpmat, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1408 CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
1409 CALL cp_dbcsr_alloc_block_from_nbl(cmatrix_db, sab_nl)
1410
1411 CALL get_kpoint_info(kpoints, mpools=mpools_kp)
1412 CALL mpools_get(mpools_kp, ao_ao_fm_pools=ao_ao_fm_pools_kp)
1413 CALL fm_pool_create_fm(ao_ao_fm_pools_kp(1)%pool, fmlocal)
1414 CALL cp_fm_get_info(fmlocal, matrix_struct=ao_ao_struct)
1415
1416 CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct)
1417 CALL cp_cfm_create(cfm_nao_nmo_work, mo_coeff%matrix_struct)
1418
1419 NULLIFY (nmo_nmo_struct)
1420 CALL cp_fm_struct_create(nmo_nmo_struct, template_fmstruct=mo_coeff%matrix_struct, &
1421 nrow_global=nmo, ncol_global=nmo)
1422 CALL cp_cfm_create(csc_cfm, nmo_nmo_struct)
1423 CALL cp_fm_struct_release(nmo_nmo_struct)
1424
1425 para_env => kpoints%blacs_env_all%para_env
1426 ALLOCATE (info(kplocal*nkp_groups, 2))
1427
1428 ALLOCATE (csmat_cur(kplocal))
1429 DO ikp = 1, kplocal
1430 CALL cp_cfm_create(csmat_cur(ikp), ao_ao_struct)
1431 END DO
1432
1433 ! Phase A: Fourier Transform S(R) -> S(k)
1434 indx = 0
1435 DO ikp = 1, kplocal
1436 DO igroup = 1, nkp_groups
1437 ik = kp_dist(1, igroup) + ikp - 1
1438 my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
1439 indx = indx + 1
1440
1441 CALL dbcsr_set(rmatrix, 0.0_dp)
1442 CALL dbcsr_set(cmatrix_db, 0.0_dp)
1443 CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix_db, rsmat=matrix_s_kp, &
1444 ispin=1, xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1445 CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1446 CALL copy_dbcsr_to_fm(tmpmat, scf_env%scf_work1(1))
1447 CALL dbcsr_desymmetrize(cmatrix_db, tmpmat)
1448 CALL copy_dbcsr_to_fm(tmpmat, scf_env%scf_work1(2))
1449
1450 IF (my_kpgrp) THEN
1451 CALL cp_fm_start_copy_general(scf_env%scf_work1(1), fmlocal, para_env, info(indx, 1))
1452 CALL cp_fm_start_copy_general(scf_env%scf_work1(2), fmlocal, para_env, info(indx, 2))
1453 ELSE
1454 CALL cp_fm_start_copy_general(scf_env%scf_work1(1), fmdummy, para_env, info(indx, 1))
1455 CALL cp_fm_start_copy_general(scf_env%scf_work1(2), fmdummy, para_env, info(indx, 2))
1456 END IF
1457 END DO
1458 END DO
1459
1460 ! Finish Communication
1461 indx = 0
1462 DO ikp = 1, kplocal
1463 DO igroup = 1, nkp_groups
1464 ik = kp_dist(1, igroup) + ikp - 1
1465 my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
1466 indx = indx + 1
1467 IF (my_kpgrp) THEN
1468 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
1469 CALL cp_cfm_scale_and_add_fm(z_zero, csmat_cur(ikp), z_one, fmlocal)
1470 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
1471 CALL cp_cfm_scale_and_add_fm(z_one, csmat_cur(ikp), gaussi, fmlocal)
1472 END IF
1473 END DO
1474 END DO
1475
1476 DO indx = 1, kplocal*nkp_groups
1477 CALL cp_fm_cleanup_copy_general(info(indx, 1))
1478 CALL cp_fm_cleanup_copy_general(info(indx, 2))
1479 END DO
1480
1481 ! Phase B: bring the WFN from its saved/internal PBC image convention to
1482 ! the current convention, then orthogonalize it with respect to S(k).
1483 ALLOCATE (eigenvalues(nmo))
1484 eval_thresh = 1.0e-12_dp
1485
1486 DO ikp = 1, kplocal
1487 kp => kpoints%kp_env(ikp)%kpoint_env
1488 DO ispin = 1, nspin
1489 CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
1490 CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1491 IF (reload_snapshot_wf) THEN
1492 CALL cp_fm_to_fm(t1_state%wf_kp(ikp, 1, ispin), rmos)
1493 CALL cp_fm_to_fm(t1_state%wf_kp(ikp, 2, ispin), imos)
1494 END IF
1495 IF (use_pbc_phase_ref) THEN
1496 ik = kp_range(1) + ikp - 1
1497 CALL wfi_apply_kp_pbc_phase_fm(rmos, imos, pbc_shift_cur - pbc_shift_src, &
1498 xkp(1:3, ik), matrix_s_kp(1, 1)%matrix)
1499 END IF
1500 CALL cp_fm_to_cfm(rmos, imos, cmos_new)
1501
1502 CALL cp_cfm_gemm('N', 'N', nao, nmo, nao, z_one, &
1503 csmat_cur(ikp), cmos_new, z_zero, cfm_nao_nmo_work)
1504 CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, &
1505 cmos_new, cfm_nao_nmo_work, z_zero, csc_cfm)
1506
1507 CALL cp_cfm_cholesky_decompose(csc_cfm, info_out=chol_info)
1508 IF (chol_info == 0) THEN
1509 CALL cp_cfm_triangular_multiply(csc_cfm, cmos_new, side='R', invert_tr=.true., uplo_tr='U')
1510 ELSE
1511 CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, cmos_new, cfm_nao_nmo_work, z_zero, csc_cfm)
1512 CALL cp_cfm_create(cfm_evecs, csc_cfm%matrix_struct)
1513 CALL cp_cfm_create(cfm_mhalf, csc_cfm%matrix_struct)
1514 CALL cp_cfm_heevd(csc_cfm, cfm_evecs, eigenvalues)
1515 CALL cp_cfm_to_cfm(cfm_evecs, cfm_mhalf)
1516 ALLOCATE (col_scaling(nmo))
1517 DO j = 1, nmo
1518 IF (eigenvalues(j) > eval_thresh) THEN
1519 col_scaling(j) = cmplx(1.0_dp/sqrt(eigenvalues(j)), 0.0_dp, kind=dp)
1520 ELSE
1521 col_scaling(j) = z_zero
1522 END IF
1523 END DO
1524 CALL cp_cfm_column_scale(cfm_mhalf, col_scaling)
1525 DEALLOCATE (col_scaling)
1526 CALL cp_cfm_gemm('N', 'C', nmo, nmo, nmo, z_one, cfm_mhalf, cfm_evecs, z_zero, csc_cfm)
1527 CALL cp_cfm_gemm('N', 'N', nao, nmo, nmo, z_one, cmos_new, csc_cfm, z_zero, cfm_nao_nmo_work)
1528 CALL cp_cfm_to_cfm(cfm_nao_nmo_work, cmos_new)
1529 CALL cp_cfm_release(cfm_evecs)
1530 CALL cp_cfm_release(cfm_mhalf)
1531 END IF
1532 CALL cp_cfm_to_fm(cmos_new, rmos, imos)
1533 END DO
1534 END DO
1535 DEALLOCATE (eigenvalues)
1536
1537 ! Phase C: Rebuild Density Matrix P(R)
1538 CALL kpoint_set_mo_occupation(kpoints, scf_control%smear)
1539 CALL kpoint_density_matrices(kpoints)
1540 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1541 CALL kpoint_density_transform(kpoints, rho_ao_kp, .false., &
1542 matrix_s_kp(1, 1)%matrix, sab_nl, scf_env%scf_work1)
1543 CALL qs_rho_update_rho(rho, qs_env=qs_env)
1544 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
1545
1546 ! Cleanup
1547 DO ikp = 1, kplocal
1548 CALL cp_cfm_release(csmat_cur(ikp))
1549 END DO
1550 DEALLOCATE (csmat_cur)
1551 DEALLOCATE (info)
1552 CALL cp_cfm_release(cmos_new)
1553 CALL cp_cfm_release(cfm_nao_nmo_work)
1554 CALL cp_cfm_release(csc_cfm)
1555 CALL fm_pool_give_back_fm(ao_ao_fm_pools_kp(1)%pool, fmlocal)
1556 CALL dbcsr_deallocate_matrix(rmatrix)
1557 CALL dbcsr_deallocate_matrix(cmatrix_db)
1558 CALL dbcsr_deallocate_matrix(tmpmat)
1559 IF (ALLOCATED(pbc_shift_cur)) DEALLOCATE (pbc_shift_cur)
1560 IF (ALLOCATED(pbc_shift_src)) DEALLOCATE (pbc_shift_src)
1561
1562 CALL timestop(handle)
1563 END SUBROUTINE wfi_use_prev_wf_kp
1564
1565! **************************************************************************************************
1566!> \brief Stores the internal PBC image shift used for k-point neighbor-list construction.
1567!> shift = scaled(pbc(r))-scaled(r), i.e. the integer image displacement caused by pbc().
1568!> \param snapshot ...
1569!> \param cell ...
1570!> \param particle_set ...
1571! **************************************************************************************************
1572 SUBROUTINE wfi_store_kp_pbc_shift(snapshot, cell, particle_set)
1573 TYPE(qs_wf_snapshot_type), POINTER :: snapshot
1574 TYPE(cell_type), POINTER :: cell
1575 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1576
1577 INTEGER :: iatom, natom
1578 REAL(kind=dp), DIMENSION(3) :: frac_pbc, frac_raw, r_pbc
1579
1580 cpassert(ASSOCIATED(snapshot))
1581 cpassert(ASSOCIATED(cell))
1582 cpassert(ASSOCIATED(particle_set))
1583
1584 natom = SIZE(particle_set)
1585 IF (ASSOCIATED(snapshot%kp_pbc_shift)) THEN
1586 DEALLOCATE (snapshot%kp_pbc_shift)
1587 END IF
1588 ALLOCATE (snapshot%kp_pbc_shift(3, natom))
1589 DO iatom = 1, natom
1590 r_pbc(1:3) = pbc(particle_set(iatom)%r(1:3), cell)
1591 CALL real_to_scaled(frac_raw, particle_set(iatom)%r(1:3), cell)
1592 CALL real_to_scaled(frac_pbc, r_pbc(1:3), cell)
1593 snapshot%kp_pbc_shift(1:3, iatom) = nint(frac_pbc(1:3) - frac_raw(1:3))
1594 END DO
1595 END SUBROUTINE wfi_store_kp_pbc_shift
1596
1597! **************************************************************************************************
1598!> \brief Computes the current internal PBC image shift used by pbc().
1599!> \param qs_env ...
1600!> \param pbc_shift ...
1601! **************************************************************************************************
1602 SUBROUTINE wfi_compute_kp_pbc_shift(qs_env, pbc_shift)
1603 TYPE(qs_environment_type), POINTER :: qs_env
1604 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: pbc_shift
1605
1606 INTEGER :: iatom, natom
1607 REAL(kind=dp), DIMENSION(3) :: frac_pbc, frac_raw, r_pbc
1608 TYPE(cell_type), POINTER :: cell
1609 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1610
1611 NULLIFY (cell, particle_set)
1612 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1613 cpassert(ASSOCIATED(cell))
1614 cpassert(ASSOCIATED(particle_set))
1615
1616 natom = SIZE(particle_set)
1617 ALLOCATE (pbc_shift(3, natom))
1618 DO iatom = 1, natom
1619 r_pbc(1:3) = pbc(particle_set(iatom)%r(1:3), cell)
1620 CALL real_to_scaled(frac_raw, particle_set(iatom)%r(1:3), cell)
1621 CALL real_to_scaled(frac_pbc, r_pbc(1:3), cell)
1622 pbc_shift(1:3, iatom) = nint(frac_pbc(1:3) - frac_raw(1:3))
1623 END DO
1624 END SUBROUTINE wfi_compute_kp_pbc_shift
1625
1626! **************************************************************************************************
1627!> \brief Applies the atom-wise Bloch phase associated with a change of the internal
1628!> k-point PBC image convention to real/imaginary MO coefficient matrices.
1629!> \param rmos real part of the MO coefficients
1630!> \param imos imaginary part of the MO coefficients
1631!> \param pbc_shift_delta target shift minus source shift for each atom
1632!> \param xk fractional k-point coordinates
1633!> \param matrix_template AO block structure used to map rows to atoms
1634! **************************************************************************************************
1635 SUBROUTINE wfi_apply_kp_pbc_phase_fm(rmos, imos, pbc_shift_delta, xk, matrix_template)
1636 TYPE(cp_fm_type), POINTER :: rmos, imos
1637 INTEGER, DIMENSION(:, :), INTENT(IN) :: pbc_shift_delta
1638 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: xk
1639 TYPE(dbcsr_type), POINTER :: matrix_template
1640
1641 INTEGER :: iatom, icol, irow, natom, nmo, nrow, &
1642 row_start
1643 INTEGER, DIMENSION(:), POINTER :: row_blk_size
1644 REAL(kind=dp) :: ci, cr, i_old, r_old, theta
1645 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: iblock, rblock
1646
1647 cpassert(ASSOCIATED(rmos))
1648 cpassert(ASSOCIATED(imos))
1649 cpassert(ASSOCIATED(matrix_template))
1650
1651 natom = SIZE(pbc_shift_delta, 2)
1652 CALL cp_fm_get_info(rmos, ncol_global=nmo)
1653 NULLIFY (row_blk_size)
1654 CALL dbcsr_get_info(matrix_template, row_blk_size=row_blk_size)
1655 cpassert(SIZE(row_blk_size) >= natom)
1656
1657 row_start = 1
1658 DO iatom = 1, natom
1659 nrow = row_blk_size(iatom)
1660 IF (any(pbc_shift_delta(1:3, iatom) /= 0)) THEN
1661 theta = twopi*sum(xk(1:3)*real(pbc_shift_delta(1:3, iatom), kind=dp))
1662 cr = cos(theta)
1663 ci = sin(theta)
1664 ALLOCATE (rblock(nrow, nmo), iblock(nrow, nmo))
1665 CALL cp_fm_get_submatrix(rmos, rblock, row_start, 1, nrow, nmo)
1666 CALL cp_fm_get_submatrix(imos, iblock, row_start, 1, nrow, nmo)
1667 DO icol = 1, nmo
1668 DO irow = 1, nrow
1669 r_old = rblock(irow, icol)
1670 i_old = iblock(irow, icol)
1671 rblock(irow, icol) = cr*r_old - ci*i_old
1672 iblock(irow, icol) = ci*r_old + cr*i_old
1673 END DO
1674 END DO
1675 CALL cp_fm_set_submatrix(rmos, rblock, row_start, 1, nrow, nmo)
1676 CALL cp_fm_set_submatrix(imos, iblock, row_start, 1, nrow, nmo)
1677 DEALLOCATE (rblock, iblock)
1678 END IF
1679 row_start = row_start + nrow
1680 END DO
1681 END SUBROUTINE wfi_apply_kp_pbc_phase_fm
1682
1683! **************************************************************************************************
1684!> \brief Applies the atom-wise Bloch phase associated with a change of the internal
1685!> k-point PBC image convention to a complex MO coefficient matrix.
1686!> \param cmos complex MO coefficients
1687!> \param pbc_shift_delta target shift minus source shift for each atom
1688!> \param xk fractional k-point coordinates
1689!> \param matrix_template AO block structure used to map rows to atoms
1690! **************************************************************************************************
1691 SUBROUTINE wfi_apply_kp_pbc_phase_cfm(cmos, pbc_shift_delta, xk, matrix_template)
1692 TYPE(cp_cfm_type), INTENT(INOUT) :: cmos
1693 INTEGER, DIMENSION(:, :), INTENT(IN) :: pbc_shift_delta
1694 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: xk
1695 TYPE(dbcsr_type), POINTER :: matrix_template
1696
1697 COMPLEX(KIND=dp) :: phase
1698 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: zblock
1699 INTEGER :: iatom, natom, nmo, nrow, row_start
1700 INTEGER, DIMENSION(:), POINTER :: row_blk_size
1701 REAL(kind=dp) :: theta
1702
1703 cpassert(ASSOCIATED(matrix_template))
1704
1705 natom = SIZE(pbc_shift_delta, 2)
1706 CALL cp_cfm_get_info(cmos, ncol_global=nmo)
1707 NULLIFY (row_blk_size)
1708 CALL dbcsr_get_info(matrix_template, row_blk_size=row_blk_size)
1709 cpassert(SIZE(row_blk_size) >= natom)
1710
1711 row_start = 1
1712 DO iatom = 1, natom
1713 nrow = row_blk_size(iatom)
1714 IF (any(pbc_shift_delta(1:3, iatom) /= 0)) THEN
1715 theta = twopi*sum(xk(1:3)*real(pbc_shift_delta(1:3, iatom), kind=dp))
1716 phase = cmplx(cos(theta), sin(theta), kind=dp)
1717 ALLOCATE (zblock(nrow, nmo))
1718 CALL cp_cfm_get_submatrix(cmos, zblock, row_start, 1, nrow, nmo)
1719 zblock = phase*zblock
1720 CALL cp_cfm_set_submatrix(cmos, zblock, row_start, 1, nrow, nmo)
1721 DEALLOCATE (zblock)
1722 END IF
1723 row_start = row_start + nrow
1724 END DO
1725 END SUBROUTINE wfi_apply_kp_pbc_phase_cfm
1726
1727! **************************************************************************************************
1728!> \brief Performs PS/ASPC wavefunction extrapolation for k-point calculations.
1729!> Applies PS/ASPC coefficients to complex MO coefficients at each k-point,
1730!> with subspace alignment via historical overlap matrices.
1731!> Delegates final orthogonalization and density building to wfi_use_prev_wf_kp.
1732!> \param wf_history wavefunction history buffer
1733!> \param qs_env QS environment
1734!> \param nvec number of history snapshots to use
1735!> \param io_unit output unit for logging
1736!> \param print_level current print level
1737! **************************************************************************************************
1738 SUBROUTINE wfi_extrapolate_ps_aspc_kp(wf_history, qs_env, nvec, io_unit, print_level)
1739 TYPE(qs_wf_history_type), POINTER :: wf_history
1740 TYPE(qs_environment_type), POINTER :: qs_env
1741 INTEGER, INTENT(IN) :: nvec, io_unit, print_level
1742
1743 CHARACTER(len=*), PARAMETER :: routinen = 'wfi_extrapolate_ps_aspc_kp'
1744
1745 INTEGER :: handle, i, ik, ikp, ispin, kplocal, &
1746 method_nr, nao, nmo, nspin
1747 INTEGER, DIMENSION(2) :: kp_range
1748 LOGICAL :: use_real_wfn
1749 REAL(kind=dp) :: alpha_coeff
1750 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
1751 TYPE(cp_cfm_type) :: cfm_nao_nmo_work, cmos_1, cmos_i, &
1752 cmos_new, csc_cfm
1753 TYPE(cp_fm_struct_type), POINTER :: nmo_nmo_struct
1754 TYPE(cp_fm_type), POINTER :: imos, mo_coeff, rmos
1755 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp
1756 TYPE(kpoint_env_type), POINTER :: kp
1757 TYPE(kpoint_type), POINTER :: kpoints
1758 TYPE(qs_wf_snapshot_type), POINTER :: t0_state, t1_state
1759
1760 method_nr = wf_history%interpolation_method_nr
1761
1762 CALL timeset(routinen, handle)
1763 NULLIFY (kpoints, kp, mo_coeff, rmos, imos, t0_state, t1_state, nmo_nmo_struct, &
1764 matrix_s_kp, xkp)
1765
1766 CALL get_qs_env(qs_env, kpoints=kpoints, matrix_s_kp=matrix_s_kp)
1767 CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn, kp_range=kp_range, xkp=xkp)
1768 kplocal = kp_range(2) - kp_range(1) + 1
1769
1770 IF (use_real_wfn) THEN
1771 IF (method_nr == wfi_aspc_nr) THEN
1772 CALL cp_warn(__location__, "ASPC with k-points requires complex wavefunctions; "// &
1773 "falling back to USE_PREV_WF.")
1774 ELSE
1775 CALL cp_warn(__location__, "PS with k-points requires complex wavefunctions; "// &
1776 "falling back to USE_PREV_WF.")
1777 END IF
1778 CALL wfi_use_prev_wf_kp(qs_env, io_unit, print_level)
1779 CALL timestop(handle)
1780 RETURN
1781 END IF
1782
1783 kp => kpoints%kp_env(1)%kpoint_env
1784 nspin = SIZE(kp%mos, 2)
1785 CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo, mo_coeff=mo_coeff)
1786
1787 IF (method_nr == wfi_aspc_nr) THEN
1788 IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1789 WRITE (unit=io_unit, fmt="(/,T2,A,/,T3,A,I0)") &
1790 "Parameters for the always stable predictor-corrector (ASPC) method:", &
1791 "ASPC order: ", max(nvec - 2, 0)
1792 END IF
1793 END IF
1794
1795 IF (method_nr == wfi_aspc_nr) THEN
1796 CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct, set_zero=.true.)
1797 CALL cp_cfm_create(cmos_1, mo_coeff%matrix_struct, set_zero=.true.)
1798 CALL cp_cfm_create(cmos_i, mo_coeff%matrix_struct, set_zero=.true.)
1799 CALL cp_cfm_create(cfm_nao_nmo_work, mo_coeff%matrix_struct, set_zero=.true.)
1800 ELSE
1801 CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct)
1802 CALL cp_cfm_create(cmos_1, mo_coeff%matrix_struct)
1803 CALL cp_cfm_create(cmos_i, mo_coeff%matrix_struct)
1804 CALL cp_cfm_create(cfm_nao_nmo_work, mo_coeff%matrix_struct)
1805 END IF
1806
1807 CALL cp_fm_struct_create(nmo_nmo_struct, template_fmstruct=mo_coeff%matrix_struct, &
1808 nrow_global=nmo, ncol_global=nmo)
1809 IF (method_nr == wfi_aspc_nr) THEN
1810 CALL cp_cfm_create(csc_cfm, nmo_nmo_struct, set_zero=.true.)
1811 ELSE
1812 CALL cp_cfm_create(csc_cfm, nmo_nmo_struct)
1813 END IF
1814 CALL cp_fm_struct_release(nmo_nmo_struct)
1815
1816 ! Phase 1: Initialize C_new(k) = B(1) * C_1(k)
1817 t1_state => wfi_get_snapshot(wf_history, wf_index=1)
1818 IF (method_nr == wfi_aspc_nr) THEN
1819 alpha_coeff = real(4*nvec - 2, kind=dp)/real(nvec + 1, kind=dp)
1820 IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1821 WRITE (unit=io_unit, fmt="(T3,A2,I0,A4,F10.6)") "B(", 1, ") = ", alpha_coeff
1822 END IF
1823 ELSE
1824 alpha_coeff = nvec
1825 END IF
1826
1827 DO ikp = 1, kplocal
1828 kp => kpoints%kp_env(ikp)%kpoint_env
1829 DO ispin = 1, nspin
1830 CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
1831 CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1832 CALL cp_fm_to_fm(t1_state%wf_kp(ikp, 1, ispin), rmos)
1833 CALL cp_fm_to_fm(t1_state%wf_kp(ikp, 2, ispin), imos)
1834 CALL cp_fm_scale(alpha_coeff, rmos)
1835 CALL cp_fm_scale(alpha_coeff, imos)
1836 END DO
1837 END DO
1838
1839 ! Phase 2: Accumulate historical snapshots C_new += B(i) * C_proj(k)
1840 DO i = 2, nvec
1841 t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1842 IF (method_nr == wfi_aspc_nr) THEN
1843 alpha_coeff = (-1.0_dp)**(i + 1)*real(i, kind=dp)* &
1844 binomial(2*nvec, nvec - i)/binomial(2*nvec - 2, nvec - 1)
1845 IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1846 WRITE (unit=io_unit, fmt="(T3,A2,I0,A4,F10.6)") "B(", i, ") = ", alpha_coeff
1847 END IF
1848 ELSE
1849 alpha_coeff = -1.0_dp*alpha_coeff*real(nvec - i + 1, dp)/real(i, dp)
1850 END IF
1851
1852 DO ikp = 1, kplocal
1853 kp => kpoints%kp_env(ikp)%kpoint_env
1854 DO ispin = 1, nspin
1855 ik = kp_range(1) + ikp - 1
1856 CALL cp_fm_to_cfm(t1_state%wf_kp(ikp, 1, ispin), t1_state%wf_kp(ikp, 2, ispin), cmos_1)
1857 CALL cp_fm_to_cfm(t0_state%wf_kp(ikp, 1, ispin), t0_state%wf_kp(ikp, 2, ispin), cmos_i)
1858
1859 ! Express the reference snapshot in the image convention of snapshot i,
1860 ! because the historical overlap below belongs to snapshot i.
1861 CALL wfi_apply_kp_pbc_phase_cfm(cmos_1, t0_state%kp_pbc_shift - t1_state%kp_pbc_shift, &
1862 xkp(1:3, ik), matrix_s_kp(1, 1)%matrix)
1863
1864 ! Subspace projection: C_proj = C_i * (C_i^dag S_i C_1)
1865 CALL cp_cfm_gemm('N', 'N', nao, nmo, nao, z_one, &
1866 t0_state%overlap_cfm_kp(ikp), cmos_1, z_zero, cfm_nao_nmo_work)
1867 CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, &
1868 cmos_i, cfm_nao_nmo_work, z_zero, csc_cfm)
1869 CALL cp_cfm_gemm('N', 'N', nao, nmo, nmo, z_one, &
1870 cmos_i, csc_cfm, z_zero, cfm_nao_nmo_work)
1871
1872 ! Convert the projected contribution from snapshot i to the reference
1873 ! image convention of snapshot 1. The final conversion to the current
1874 ! convention is centralized in wfi_use_prev_wf_kp.
1875 CALL wfi_apply_kp_pbc_phase_cfm(cfm_nao_nmo_work, t1_state%kp_pbc_shift - t0_state%kp_pbc_shift, &
1876 xkp(1:3, ik), matrix_s_kp(1, 1)%matrix)
1877
1878 CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
1879 CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1880 CALL cp_fm_to_cfm(rmos, imos, cmos_new)
1881 CALL cp_cfm_scale_and_add(z_one, cmos_new, cmplx(alpha_coeff, 0.0_dp, kind=dp), cfm_nao_nmo_work)
1882 CALL cp_cfm_to_fm(cmos_new, rmos, imos)
1883 END DO
1884 END DO
1885 END DO
1886
1887 CALL cp_cfm_release(cmos_new)
1888 CALL cp_cfm_release(cmos_1)
1889 CALL cp_cfm_release(cmos_i)
1890 CALL cp_cfm_release(cfm_nao_nmo_work)
1891 CALL cp_cfm_release(csc_cfm)
1892
1893 ! Phase 3: Convert the extrapolated WFN from the reference snapshot image
1894 ! convention to the current k-point PBC convention, then reorthogonalize and
1895 ! rebuild the density. Keep the actual phase handling centralized in
1896 ! wfi_use_prev_wf_kp so that USE_PREV_WF and ASPC/PS share the same path.
1897 CALL wfi_use_prev_wf_kp(qs_env, 0, print_level, pbc_shift_ref=t1_state%kp_pbc_shift, &
1898 load_snapshot_wf=.false.)
1899
1900 CALL timestop(handle)
1901
1902 END SUBROUTINE wfi_extrapolate_ps_aspc_kp
1903
1904! **************************************************************************************************
1905!> \brief GEXT_PROJ/GEXT_PROJ_QTR wavefunction extrapolation for complex k-points.
1906!> This follows the existing ASPC/PS k-point projection path, but uses
1907!> the GEXT-fitted coefficients.
1908!> \param wf_history wavefunction history buffer
1909!> \param qs_env The QS environment
1910!> \param nvec number of previous wavefunctions
1911!> \param io_unit output unit
1912!> \param print_level current print level
1913! **************************************************************************************************
1914 SUBROUTINE wfi_extrapolate_gext_proj_kp(wf_history, qs_env, nvec, io_unit, print_level)
1915 TYPE(qs_wf_history_type), POINTER :: wf_history
1916 TYPE(qs_environment_type), POINTER :: qs_env
1917 INTEGER, INTENT(IN) :: nvec, io_unit, print_level
1918
1919 CHARACTER(len=*), PARAMETER :: routinen = 'wfi_extrapolate_gext_proj_kp'
1920
1921 INTEGER :: handle, i, igroup, ik, ikp, indx, ispin, &
1922 kplocal, method_nr, nao, nkp, &
1923 nkp_groups, nmo, nspin
1924 INTEGER, DIMENSION(2) :: kp_range
1925 INTEGER, DIMENSION(:, :), POINTER :: kp_dist
1926 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1927 LOGICAL :: my_kpgrp, use_real_wfn
1928 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coeffs, weight_kp
1929 REAL(kind=dp), DIMENSION(:), POINTER :: wkp
1930 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
1931 TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
1932 TYPE(cp_cfm_type) :: cfm_nao_nmo_work, cmos_1, cmos_i, &
1933 cmos_new, csc_cfm
1934 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: csmat_cur
1935 TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools_kp
1936 TYPE(cp_fm_struct_type), POINTER :: ao_ao_struct, nmo_nmo_struct
1937 TYPE(cp_fm_type) :: fmdummy, fmlocal
1938 TYPE(cp_fm_type), POINTER :: imos, mo_coeff, rmos
1939 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp
1940 TYPE(dbcsr_type), POINTER :: cmatrix_db, rmatrix, tmpmat
1941 TYPE(kpoint_env_type), POINTER :: kp
1942 TYPE(kpoint_type), POINTER :: kpoints
1943 TYPE(mp_para_env_type), POINTER :: para_env, para_env_inter_kp
1944 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1945 POINTER :: sab_nl
1946 TYPE(qs_matrix_pools_type), POINTER :: mpools_kp
1947 TYPE(qs_scf_env_type), POINTER :: scf_env
1948 TYPE(qs_wf_snapshot_type), POINTER :: t0_state, t1_state
1949
1950 method_nr = wf_history%interpolation_method_nr
1951
1952 CALL timeset(routinen, handle)
1953 NULLIFY (ao_ao_struct, cell_to_index, cmatrix_db, imos, kp, kpoints, matrix_s_kp, &
1954 mo_coeff, mpools_kp, para_env, para_env_inter_kp, rmatrix, rmos, sab_nl, &
1955 scf_env, t0_state, t1_state, tmpmat, xkp, kp_dist, wkp, nmo_nmo_struct, &
1956 ao_ao_fm_pools_kp)
1957
1958 CALL get_qs_env(qs_env, kpoints=kpoints, matrix_s_kp=matrix_s_kp, scf_env=scf_env)
1959 CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn, kp_range=kp_range, &
1960 nkp=nkp, xkp=xkp, wkp=wkp, nkp_groups=nkp_groups, &
1961 kp_dist=kp_dist, cell_to_index=cell_to_index, sab_nl=sab_nl, &
1962 mpools=mpools_kp, para_env_inter_kp=para_env_inter_kp)
1963 kplocal = kp_range(2) - kp_range(1) + 1
1964
1965 IF (use_real_wfn) THEN
1966 CALL cp_warn(__location__, "GExt with k-points requires complex wavefunctions; "// &
1967 "falling back to USE_PREV_WF.")
1968 CALL wfi_use_prev_wf_kp(qs_env, io_unit, print_level)
1969 CALL timestop(handle)
1970 RETURN
1971 END IF
1972
1973 IF (nvec >= wf_history%memory_depth) THEN
1974 IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. &
1975 (qs_env%scf_control%eps_scf_hist /= 0)) THEN
1976 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1977 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1978 qs_env%scf_control%outer_scf%have_scf = .false.
1979 ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
1980 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1981 qs_env%scf_control%outer_scf%have_scf = .false.
1982 ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
1983 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1984 END IF
1985 END IF
1986
1987 kp => kpoints%kp_env(1)%kpoint_env
1988 nspin = SIZE(kp%mos, 2)
1989 CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo, mo_coeff=mo_coeff)
1990
1991 ! Build the current S(k), using the same Fourier-transform pattern as
1992 ! wfi_use_prev_wf_kp. This is only needed for the GEXT coefficient fitting.
1993 ALLOCATE (rmatrix, cmatrix_db, tmpmat)
1994 CALL dbcsr_create(rmatrix, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
1995 CALL dbcsr_create(cmatrix_db, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
1996 CALL dbcsr_create(tmpmat, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1997 CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
1998 CALL cp_dbcsr_alloc_block_from_nbl(cmatrix_db, sab_nl)
1999
2000 CALL mpools_get(mpools_kp, ao_ao_fm_pools=ao_ao_fm_pools_kp)
2001 CALL fm_pool_create_fm(ao_ao_fm_pools_kp(1)%pool, fmlocal)
2002 CALL cp_fm_get_info(fmlocal, matrix_struct=ao_ao_struct)
2003
2004 para_env => kpoints%blacs_env_all%para_env
2005 ALLOCATE (info(kplocal*nkp_groups, 2))
2006 ALLOCATE (csmat_cur(kplocal), weight_kp(kplocal))
2007 DO ikp = 1, kplocal
2008 CALL cp_cfm_create(csmat_cur(ikp), ao_ao_struct)
2009 weight_kp(ikp) = wkp(kp_range(1) + ikp - 1)
2010 END DO
2011
2012 indx = 0
2013 DO ikp = 1, kplocal
2014 DO igroup = 1, nkp_groups
2015 ik = kp_dist(1, igroup) + ikp - 1
2016 my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
2017 indx = indx + 1
2018
2019 CALL dbcsr_set(rmatrix, 0.0_dp)
2020 CALL dbcsr_set(cmatrix_db, 0.0_dp)
2021 CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix_db, rsmat=matrix_s_kp, &
2022 ispin=1, xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
2023 CALL dbcsr_desymmetrize(rmatrix, tmpmat)
2024 CALL copy_dbcsr_to_fm(tmpmat, scf_env%scf_work1(1))
2025 CALL dbcsr_desymmetrize(cmatrix_db, tmpmat)
2026 CALL copy_dbcsr_to_fm(tmpmat, scf_env%scf_work1(2))
2027
2028 IF (my_kpgrp) THEN
2029 CALL cp_fm_start_copy_general(scf_env%scf_work1(1), fmlocal, para_env, info(indx, 1))
2030 CALL cp_fm_start_copy_general(scf_env%scf_work1(2), fmlocal, para_env, info(indx, 2))
2031 ELSE
2032 CALL cp_fm_start_copy_general(scf_env%scf_work1(1), fmdummy, para_env, info(indx, 1))
2033 CALL cp_fm_start_copy_general(scf_env%scf_work1(2), fmdummy, para_env, info(indx, 2))
2034 END IF
2035 END DO
2036 END DO
2037
2038 indx = 0
2039 DO ikp = 1, kplocal
2040 DO igroup = 1, nkp_groups
2041 ik = kp_dist(1, igroup) + ikp - 1
2042 my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
2043 indx = indx + 1
2044 IF (my_kpgrp) THEN
2045 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
2046 CALL cp_cfm_scale_and_add_fm(z_zero, csmat_cur(ikp), z_one, fmlocal)
2047 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
2048 CALL cp_cfm_scale_and_add_fm(z_one, csmat_cur(ikp), gaussi, fmlocal)
2049 END IF
2050 END DO
2051 END DO
2052
2053 DO indx = 1, kplocal*nkp_groups
2054 CALL cp_fm_cleanup_copy_general(info(indx, 1))
2055 CALL cp_fm_cleanup_copy_general(info(indx, 2))
2056 END DO
2057
2058 ALLOCATE (coeffs(nvec))
2059 IF (method_nr == wfi_gext_proj_nr) THEN
2060 CALL diff_fitting(wf_history, matrix_s_kp(1, 1)%matrix, coeffs, nvec, &
2061 1e-4_dp, io_unit, print_level, current_overlap_kp=csmat_cur, &
2062 kpoint_weights=weight_kp, para_env_inter_kp=para_env_inter_kp)
2063 ELSE
2064 CALL tr_fitting(wf_history, matrix_s_kp(1, 1)%matrix, coeffs, nvec, &
2065 1e-4_dp, io_unit, print_level, current_overlap_kp=csmat_cur, &
2066 kpoint_weights=weight_kp, para_env_inter_kp=para_env_inter_kp)
2067 END IF
2068
2069 ! Accumulate the extrapolated WFN using the same projected-WFN path as ASPC/PS.
2070 CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct)
2071 CALL cp_cfm_create(cmos_1, mo_coeff%matrix_struct)
2072 CALL cp_cfm_create(cmos_i, mo_coeff%matrix_struct)
2073 CALL cp_cfm_create(cfm_nao_nmo_work, mo_coeff%matrix_struct)
2074 CALL cp_fm_struct_create(nmo_nmo_struct, template_fmstruct=mo_coeff%matrix_struct, &
2075 nrow_global=nmo, ncol_global=nmo)
2076 CALL cp_cfm_create(csc_cfm, nmo_nmo_struct)
2077 CALL cp_fm_struct_release(nmo_nmo_struct)
2078
2079 t1_state => wfi_get_snapshot(wf_history, wf_index=1)
2080 DO ikp = 1, kplocal
2081 kp => kpoints%kp_env(ikp)%kpoint_env
2082 DO ispin = 1, nspin
2083 CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
2084 CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
2085 CALL cp_fm_set_all(rmos, 0.0_dp)
2086 CALL cp_fm_set_all(imos, 0.0_dp)
2087 END DO
2088 END DO
2089
2090 DO i = 1, nvec
2091 t0_state => wfi_get_snapshot(wf_history, wf_index=i)
2092 DO ikp = 1, kplocal
2093 kp => kpoints%kp_env(ikp)%kpoint_env
2094 ik = kp_range(1) + ikp - 1
2095 DO ispin = 1, nspin
2096 CALL cp_fm_to_cfm(t1_state%wf_kp(ikp, 1, ispin), t1_state%wf_kp(ikp, 2, ispin), cmos_1)
2097 CALL cp_fm_to_cfm(t0_state%wf_kp(ikp, 1, ispin), t0_state%wf_kp(ikp, 2, ispin), cmos_i)
2098
2099 CALL wfi_apply_kp_pbc_phase_cfm(cmos_1, t0_state%kp_pbc_shift - t1_state%kp_pbc_shift, &
2100 xkp(1:3, ik), matrix_s_kp(1, 1)%matrix)
2101
2102 CALL cp_cfm_gemm('N', 'N', nao, nmo, nao, z_one, &
2103 t0_state%overlap_cfm_kp(ikp), cmos_1, z_zero, cfm_nao_nmo_work)
2104 CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, &
2105 cmos_i, cfm_nao_nmo_work, z_zero, csc_cfm)
2106 CALL cp_cfm_gemm('N', 'N', nao, nmo, nmo, z_one, &
2107 cmos_i, csc_cfm, z_zero, cfm_nao_nmo_work)
2108
2109 CALL wfi_apply_kp_pbc_phase_cfm(cfm_nao_nmo_work, t1_state%kp_pbc_shift - t0_state%kp_pbc_shift, &
2110 xkp(1:3, ik), matrix_s_kp(1, 1)%matrix)
2111
2112 CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
2113 CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
2114 CALL cp_fm_to_cfm(rmos, imos, cmos_new)
2115 CALL cp_cfm_scale_and_add(z_one, cmos_new, cmplx(coeffs(i), 0.0_dp, kind=dp), cfm_nao_nmo_work)
2116 CALL cp_cfm_to_fm(cmos_new, rmos, imos)
2117 END DO
2118 END DO
2119 END DO
2120
2121 CALL cp_cfm_release(cmos_new)
2122 CALL cp_cfm_release(cmos_1)
2123 CALL cp_cfm_release(cmos_i)
2124 CALL cp_cfm_release(cfm_nao_nmo_work)
2125 CALL cp_cfm_release(csc_cfm)
2126
2127 CALL wfi_use_prev_wf_kp(qs_env, 0, print_level, pbc_shift_ref=t1_state%kp_pbc_shift, &
2128 load_snapshot_wf=.false.)
2129
2130 DO ikp = 1, kplocal
2131 CALL cp_cfm_release(csmat_cur(ikp))
2132 END DO
2133 DEALLOCATE (csmat_cur, coeffs, weight_kp, info)
2134 CALL fm_pool_give_back_fm(ao_ao_fm_pools_kp(1)%pool, fmlocal)
2135 CALL dbcsr_deallocate_matrix(rmatrix)
2136 CALL dbcsr_deallocate_matrix(cmatrix_db)
2137 CALL dbcsr_deallocate_matrix(tmpmat)
2138
2139 CALL timestop(handle)
2140
2141 END SUBROUTINE wfi_extrapolate_gext_proj_kp
2142
2143! **************************************************************************************************
2144!> \brief Decides if scf control variables has to changed due
2145!> to using a WF extrapolation.
2146!> \param qs_env The QS environment
2147!> \param nvec ...
2148!> \par History
2149!> 11.2006 created [TdK]
2150!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch)
2151! **************************************************************************************************
2152 ELEMENTAL SUBROUTINE wfi_set_history_variables(qs_env, nvec)
2153 TYPE(qs_environment_type), INTENT(INOUT) :: qs_env
2154 INTEGER, INTENT(IN) :: nvec
2155
2156 IF (nvec >= qs_env%wf_history%memory_depth) THEN
2157 IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. (qs_env%scf_control%eps_scf_hist /= 0)) THEN
2158 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
2159 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
2160 qs_env%scf_control%outer_scf%have_scf = .false.
2161 ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
2162 qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
2163 qs_env%scf_control%outer_scf%have_scf = .false.
2164 ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
2165 qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
2166 qs_env%scf_control%outer_scf%eps_scf = qs_env%scf_control%eps_scf_hist
2167 END IF
2168 END IF
2169
2170 END SUBROUTINE wfi_set_history_variables
2171
2172! **************************************************************************************************
2173!> \brief updates the snapshot buffer, taking a new snapshot
2174!> \param wf_history the history buffer to update
2175!> \param qs_env the qs_env we get the info from
2176!> \param dt ...
2177!> \par History
2178!> 02.2003 created [fawzi]
2179!> \author fawzi
2180! **************************************************************************************************
2181 SUBROUTINE wfi_update(wf_history, qs_env, dt)
2182 TYPE(qs_wf_history_type), POINTER :: wf_history
2183 TYPE(qs_environment_type), POINTER :: qs_env
2184 REAL(kind=dp), INTENT(in) :: dt
2185
2186 cpassert(ASSOCIATED(wf_history))
2187 cpassert(wf_history%ref_count > 0)
2188 cpassert(ASSOCIATED(qs_env))
2189
2190 wf_history%snapshot_count = wf_history%snapshot_count + 1
2191 IF (wf_history%memory_depth > 0) THEN
2192 wf_history%last_state_index = modulo(wf_history%snapshot_count, &
2193 wf_history%memory_depth) + 1
2194 CALL wfs_update(snapshot=wf_history%past_states &
2195 (wf_history%last_state_index)%snapshot, wf_history=wf_history, &
2196 qs_env=qs_env, dt=dt)
2197 END IF
2198 END SUBROUTINE wfi_update
2199
2200! **************************************************************************************************
2201!> \brief reorthogonalizes the mos
2202!> \param qs_env the qs_env in which to orthogonalize
2203!> \param v_matrix the vectors to orthogonalize
2204!> \param n_col number of column of v to orthogonalize
2205!> \par History
2206!> 04.2003 created [fawzi]
2207!> \author Fawzi Mohamed
2208! **************************************************************************************************
2209 SUBROUTINE reorthogonalize_vectors(qs_env, v_matrix, n_col)
2210 TYPE(qs_environment_type), POINTER :: qs_env
2211 TYPE(cp_fm_type), INTENT(IN) :: v_matrix
2212 INTEGER, INTENT(in), OPTIONAL :: n_col
2213
2214 CHARACTER(len=*), PARAMETER :: routinen = 'reorthogonalize_vectors'
2215
2216 INTEGER :: handle, my_n_col
2217 LOGICAL :: has_unit_metric, &
2218 ortho_contains_cholesky, &
2219 smearing_is_used
2220 TYPE(cp_fm_pool_type), POINTER :: maxao_maxmo_fm_pool
2221 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
2222 TYPE(dft_control_type), POINTER :: dft_control
2223 TYPE(qs_matrix_pools_type), POINTER :: mpools
2224 TYPE(qs_scf_env_type), POINTER :: scf_env
2225 TYPE(scf_control_type), POINTER :: scf_control
2226
2227 NULLIFY (scf_env, scf_control, maxao_maxmo_fm_pool, matrix_s, mpools, dft_control)
2228 CALL timeset(routinen, handle)
2229
2230 cpassert(ASSOCIATED(qs_env))
2231
2232 CALL cp_fm_get_info(v_matrix, ncol_global=my_n_col)
2233 IF (PRESENT(n_col)) my_n_col = n_col
2234 CALL get_qs_env(qs_env, mpools=mpools, &
2235 scf_env=scf_env, &
2236 scf_control=scf_control, &
2237 matrix_s=matrix_s, &
2238 dft_control=dft_control)
2239 CALL mpools_get(mpools, maxao_maxmo_fm_pool=maxao_maxmo_fm_pool)
2240 IF (ASSOCIATED(scf_env)) THEN
2241 ortho_contains_cholesky = (scf_env%method /= ot_method_nr) .AND. &
2242 (scf_env%cholesky_method > 0) .AND. &
2243 ASSOCIATED(scf_env%ortho)
2244 ELSE
2245 ortho_contains_cholesky = .false.
2246 END IF
2247
2248 CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
2249 smearing_is_used = .false.
2250 IF (dft_control%smear) THEN
2251 smearing_is_used = .true.
2252 END IF
2253
2254 IF (has_unit_metric) THEN
2255 CALL make_basis_simple(v_matrix, my_n_col)
2256 ELSE IF (smearing_is_used) THEN
2257 CALL make_basis_lowdin(vmatrix=v_matrix, ncol=my_n_col, &
2258 matrix_s=matrix_s(1)%matrix)
2259 ELSE IF (ortho_contains_cholesky) THEN
2260 CALL make_basis_cholesky(vmatrix=v_matrix, ncol=my_n_col, &
2261 ortho=scf_env%ortho)
2262 ELSE
2263 CALL make_basis_sm(v_matrix, my_n_col, matrix_s(1)%matrix)
2264 END IF
2265 CALL timestop(handle)
2266 END SUBROUTINE reorthogonalize_vectors
2267
2268! **************************************************************************************************
2269!> \brief purges wf_history retaining only the latest snapshot
2270!> \param qs_env the qs env with the latest result, and that will contain
2271!> the purged wf_history
2272!> \par History
2273!> 05.2016 created [Nico Holmberg]
2274!> \author Nico Holmberg
2275! **************************************************************************************************
2276 SUBROUTINE wfi_purge_history(qs_env)
2277 TYPE(qs_environment_type), POINTER :: qs_env
2278
2279 CHARACTER(len=*), PARAMETER :: routinen = 'wfi_purge_history'
2280
2281 INTEGER :: handle, io_unit, print_level
2282 TYPE(cp_logger_type), POINTER :: logger
2283 TYPE(dft_control_type), POINTER :: dft_control
2284 TYPE(qs_wf_history_type), POINTER :: wf_history
2285
2286 NULLIFY (dft_control, wf_history)
2287
2288 CALL timeset(routinen, handle)
2289 logger => cp_get_default_logger()
2290 print_level = logger%iter_info%print_level
2291 io_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%SCF%PRINT%PROGRAM_RUN_INFO", &
2292 extension=".scfLog")
2293
2294 cpassert(ASSOCIATED(qs_env))
2295 cpassert(ASSOCIATED(qs_env%wf_history))
2296 cpassert(qs_env%wf_history%ref_count > 0)
2297 CALL get_qs_env(qs_env, dft_control=dft_control)
2298
2299 SELECT CASE (qs_env%wf_history%interpolation_method_nr)
2300 CASE (wfi_use_guess_method_nr, wfi_use_prev_wf_method_nr, &
2301 wfi_use_prev_p_method_nr, wfi_use_prev_rho_r_method_nr, &
2302 wfi_frozen_method_nr)
2303 ! do nothing
2304 CASE (wfi_linear_wf_method_nr, wfi_linear_p_method_nr, &
2305 wfi_linear_ps_method_nr, wfi_ps_method_nr, &
2306 wfi_aspc_nr, wfi_gext_proj_nr, wfi_gext_proj_qtr_nr)
2307 IF (qs_env%wf_history%snapshot_count >= 2) THEN
2308 IF (debug_this_module .AND. io_unit > 0) &
2309 WRITE (io_unit, fmt="(T2,A)") "QS| Purging WFN history"
2310 CALL wfi_create(wf_history, interpolation_method_nr= &
2311 dft_control%qs_control%wf_interpolation_method_nr, &
2312 extrapolation_order=dft_control%qs_control%wf_extrapolation_order, &
2313 has_unit_metric=qs_env%has_unit_metric)
2314 CALL set_qs_env(qs_env=qs_env, &
2315 wf_history=wf_history)
2316 CALL wfi_release(wf_history)
2317 CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
2318 END IF
2319 CASE DEFAULT
2320 cpabort("Unknown extrapolation method.")
2321 END SELECT
2322 CALL timestop(handle)
2323
2324 END SUBROUTINE wfi_purge_history
2325
2326! **************************************************************************************************
2327!> \brief Gives the coefficients that best approximate the new overlap
2328!> as a linear combination of the previous overlaps in the
2329!> wf_history buffer. This is done by solving
2330!> argmin_a || S_{n+1} - S_{n} - \sum_i^{nvec-1} a_i (S_{n-q+i} - S_{n}) ||^2
2331!> \param wf_history wavefunction history buffer, containing the previous overlaps
2332!> \param current_overlap current overlap in dbcsr format
2333!> \param coeffs resulting nvec coefficients
2334!> \param nvec number of previous overlaps
2335!> \param eps Tikhonov regularization
2336!> \param io_unit output unit
2337!> \param print_level print level
2338!> \param current_overlap_kp ...
2339!> \param kpoint_weights ...
2340!> \param para_env_inter_kp ...
2341!> \par History
2342!> 04.2026 created [Michele Nottoli]
2343!> \author Michele Nottoli
2344! **************************************************************************************************
2345 SUBROUTINE diff_fitting(wf_history, current_overlap, coeffs, nvec, eps, io_unit, print_level, &
2346 current_overlap_kp, kpoint_weights, para_env_inter_kp)
2347 TYPE(qs_wf_history_type), POINTER :: wf_history
2348 TYPE(dbcsr_type), INTENT(IN) :: current_overlap
2349 INTEGER, INTENT(IN) :: nvec
2350 REAL(kind=dp), INTENT(OUT) :: coeffs(nvec)
2351 REAL(kind=dp), INTENT(IN) :: eps
2352 INTEGER, INTENT(IN) :: io_unit, print_level
2353 TYPE(cp_cfm_type), DIMENSION(:), INTENT(IN), &
2354 OPTIONAL :: current_overlap_kp
2355 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: kpoint_weights
2356 TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_inter_kp
2357
2358 COMPLEX(KIND=dp) :: ztrace
2359 INTEGER :: i, icol_local, ikp, info, irow_local, j
2360 REAL(kind=dp) :: error, norm_ref, weight
2361 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: b
2362 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: a
2363 TYPE(cp_cfm_type) :: target_diff_cfm, tmp_conj_cfm, &
2364 tmp_i_cfm, tmp_j_cfm
2365 TYPE(dbcsr_type) :: target_diff, tmp_i, tmp_j, tmp_k
2366 TYPE(qs_wf_snapshot_type), POINTER :: ref_state, state
2367
2368 IF (nvec <= 0) THEN
2369 cpabort("Not enough vectors to do the fitting")
2370 ELSE IF (nvec == 1) THEN
2371 coeffs(1) = 1.0_dp
2372 RETURN
2373 END IF
2374
2375 IF (PRESENT(current_overlap_kp)) THEN
2376 ALLOCATE (a(nvec - 1, nvec - 1), b(nvec - 1))
2377 a = 0.0_dp
2378 b = 0.0_dp
2379
2380 ref_state => wfi_get_snapshot(wf_history, wf_index=1)
2381 CALL cp_cfm_create(target_diff_cfm, current_overlap_kp(1)%matrix_struct)
2382 CALL cp_cfm_create(tmp_i_cfm, current_overlap_kp(1)%matrix_struct)
2383 CALL cp_cfm_create(tmp_j_cfm, current_overlap_kp(1)%matrix_struct)
2384 CALL cp_cfm_create(tmp_conj_cfm, current_overlap_kp(1)%matrix_struct)
2385
2386 DO ikp = 1, SIZE(current_overlap_kp)
2387 weight = 1.0_dp
2388 IF (PRESENT(kpoint_weights)) weight = kpoint_weights(ikp)
2389
2390 CALL cp_cfm_to_cfm(current_overlap_kp(ikp), target_diff_cfm)
2391 CALL cp_cfm_scale_and_add(z_one, target_diff_cfm, cmplx(-1.0_dp, 0.0_dp, kind=dp), &
2392 ref_state%overlap_cfm_kp(ikp))
2393 DO i = 2, nvec
2394 state => wfi_get_snapshot(wf_history, wf_index=i)
2395 CALL cp_cfm_to_cfm(state%overlap_cfm_kp(ikp), tmp_i_cfm)
2396 CALL cp_cfm_scale_and_add(z_one, tmp_i_cfm, cmplx(-1.0_dp, 0.0_dp, kind=dp), &
2397 ref_state%overlap_cfm_kp(ikp))
2398 CALL cp_cfm_to_cfm(tmp_i_cfm, tmp_conj_cfm)
2399 DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
2400 DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
2401 tmp_conj_cfm%local_data(irow_local, icol_local) = &
2402 conjg(tmp_conj_cfm%local_data(irow_local, icol_local))
2403 END DO
2404 END DO
2405 CALL cp_cfm_trace(tmp_conj_cfm, target_diff_cfm, ztrace)
2406 b(i - 1) = b(i - 1) + weight*real(ztrace, kind=dp)
2407
2408 DO j = 2, i
2409 state => wfi_get_snapshot(wf_history, wf_index=j)
2410 CALL cp_cfm_to_cfm(state%overlap_cfm_kp(ikp), tmp_j_cfm)
2411 CALL cp_cfm_scale_and_add(z_one, tmp_j_cfm, cmplx(-1.0_dp, 0.0_dp, kind=dp), &
2412 ref_state%overlap_cfm_kp(ikp))
2413 CALL cp_cfm_trace(tmp_conj_cfm, tmp_j_cfm, ztrace)
2414 a(j - 1, i - 1) = a(j - 1, i - 1) + weight*real(ztrace, kind=dp)
2415 END DO
2416 END DO
2417 END DO
2418
2419 DO i = 2, nvec
2420 DO j = 2, i
2421 a(i - 1, j - 1) = a(j - 1, i - 1)
2422 END DO
2423 END DO
2424
2425 IF (PRESENT(para_env_inter_kp)) THEN
2426 IF (ASSOCIATED(para_env_inter_kp)) THEN
2427 CALL para_env_inter_kp%sum(a)
2428 CALL para_env_inter_kp%sum(b)
2429 END IF
2430 END IF
2431
2432 DO i = 1, nvec - 1
2433 a(i, i) = a(i, i) + eps**2
2434 END DO
2435
2436 CALL dposv('u', nvec - 1, 1, a, nvec - 1, b, nvec - 1, info)
2437 IF (info /= 0) THEN
2438 cpabort("DPOSV failed.")
2439 END IF
2440
2441 coeffs(1) = 1.0_dp - sum(b)
2442 coeffs(2:nvec) = b(:)
2443
2444 IF (print_level > low_print_level) THEN
2445 error = 0.0_dp
2446 norm_ref = 0.0_dp
2447 DO ikp = 1, SIZE(current_overlap_kp)
2448 weight = 1.0_dp
2449 IF (PRESENT(kpoint_weights)) weight = kpoint_weights(ikp)
2450 CALL cp_cfm_to_cfm(current_overlap_kp(ikp), tmp_i_cfm)
2451 DO i = 1, nvec
2452 state => wfi_get_snapshot(wf_history, wf_index=i)
2453 CALL cp_cfm_scale_and_add(z_one, tmp_i_cfm, cmplx(-coeffs(i), 0.0_dp, kind=dp), &
2454 state%overlap_cfm_kp(ikp))
2455 END DO
2456 CALL cp_cfm_to_cfm(tmp_i_cfm, tmp_conj_cfm)
2457 DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
2458 DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
2459 tmp_conj_cfm%local_data(irow_local, icol_local) = &
2460 conjg(tmp_conj_cfm%local_data(irow_local, icol_local))
2461 END DO
2462 END DO
2463 CALL cp_cfm_trace(tmp_conj_cfm, tmp_i_cfm, ztrace)
2464 error = error + weight*real(ztrace, kind=dp)
2465 CALL cp_cfm_to_cfm(ref_state%overlap_cfm_kp(ikp), tmp_conj_cfm)
2466 DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
2467 DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
2468 tmp_conj_cfm%local_data(irow_local, icol_local) = &
2469 conjg(tmp_conj_cfm%local_data(irow_local, icol_local))
2470 END DO
2471 END DO
2472 CALL cp_cfm_trace(tmp_conj_cfm, ref_state%overlap_cfm_kp(ikp), ztrace)
2473 norm_ref = norm_ref + weight*real(ztrace, kind=dp)
2474 END DO
2475 IF (PRESENT(para_env_inter_kp)) THEN
2476 IF (ASSOCIATED(para_env_inter_kp)) THEN
2477 CALL para_env_inter_kp%sum(error)
2478 CALL para_env_inter_kp%sum(norm_ref)
2479 END IF
2480 END IF
2481 IF (io_unit > 0) THEN
2482 WRITE (unit=io_unit, fmt="(/,T2,A,F20.10)") "GEXT overlap fitting error:", &
2483 sqrt(error/max(norm_ref, tiny(1.0_dp)))
2484 END IF
2485 END IF
2486
2487 CALL cp_cfm_release(target_diff_cfm)
2488 CALL cp_cfm_release(tmp_i_cfm)
2489 CALL cp_cfm_release(tmp_j_cfm)
2490 CALL cp_cfm_release(tmp_conj_cfm)
2491 DEALLOCATE (a, b)
2492 RETURN
2493 END IF
2494
2495 ALLOCATE (a(nvec - 1, nvec - 1), b(nvec - 1))
2496
2497 ! get the reference for the difference fitting
2498 ref_state => wfi_get_snapshot(wf_history, wf_index=1)
2499
2500 ! assemble the target difference
2501 CALL dbcsr_copy(target_diff, current_overlap)
2502 CALL dbcsr_add(target_diff, ref_state%overlap, 1.0_dp, -1.0_dp)
2503
2504 ! allocate tmp_k
2505 CALL dbcsr_copy(tmp_k, current_overlap)
2506
2507 ! assemble the matrix A and the RHS b
2508 DO i = 2, nvec
2509 state => wfi_get_snapshot(wf_history, wf_index=i)
2510 CALL dbcsr_copy(tmp_i, state%overlap)
2511 CALL dbcsr_add(tmp_i, ref_state%overlap, 1.0_dp, -1.0_dp)
2512 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_i, target_diff, 0.0_dp, tmp_k)
2513 CALL dbcsr_trace(tmp_k, b(i - 1))
2514
2515 DO j = 2, i
2516 state => wfi_get_snapshot(wf_history, wf_index=j)
2517 CALL dbcsr_copy(tmp_j, state%overlap)
2518 CALL dbcsr_add(tmp_j, ref_state%overlap, 1.0_dp, -1.0_dp)
2519 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_i, tmp_j, 0.0_dp, tmp_k)
2520 CALL dbcsr_trace(tmp_k, a(j - 1, i - 1))
2521 a(i - 1, j - 1) = a(j - 1, i - 1)
2522 END DO
2523 END DO
2524
2525 ! add the Tikhonov regularization
2526 DO i = 1, nvec - 1
2527 a(i, i) = a(i, i) + eps**2
2528 END DO
2529
2530 ! solve the linear system
2531 CALL dposv('u', nvec - 1, 1, a, nvec - 1, b, nvec - 1, info)
2532 IF (info /= 0) THEN
2533 cpabort("DPOSV failed.")
2534 END IF
2535
2536 ! set the coefficient for the reference snapshot
2537 coeffs(1) = 1.0_dp - sum(b)
2538 coeffs(2:nvec) = b(:)
2539
2540 ! as a consistency check, print how well the current overlap
2541 ! is approximated by the linear combination of previous overlaps
2542 IF (print_level > low_print_level) THEN
2543 CALL dbcsr_copy(tmp_i, current_overlap)
2544 DO i = 1, nvec
2545 state => wfi_get_snapshot(wf_history, wf_index=i)
2546 CALL dbcsr_add(tmp_i, state%overlap, 1.0_dp, -coeffs(i))
2547 END DO
2548 error = dbcsr_frobenius_norm(tmp_i)/dbcsr_frobenius_norm(state%overlap)
2549 IF (io_unit > 0) THEN
2550 WRITE (unit=io_unit, fmt="(/,T2,A,F20.10)") "GEXT overlap fitting error:", error
2551 END IF
2552 END IF
2553
2554 ! free the memory
2555 CALL dbcsr_release(tmp_i)
2556 CALL dbcsr_release(tmp_j)
2557 CALL dbcsr_release(tmp_k)
2558 CALL dbcsr_release(target_diff)
2559 DEALLOCATE (a, b)
2560
2561 END SUBROUTINE diff_fitting
2562
2563! **************************************************************************************************
2564!> \brief Gives the coefficients that best approximate the new overlap
2565!> as a time reversible linear combination of the previous overlaps in the
2566!> wf_history buffer. This is done by solving
2567!> argmin_a || S_{n+1} + S_{n+1-nvec}
2568!> - \sum_{i=1}^q a_i (S_{n+1-nvec+i} + S_{n+1-i}) ||^2
2569!> with q = nvec/2 if nvec is even, or q = (nvec-1)/2 if odd.
2570!> \param wf_history wavefunction history buffer, containing the previous overlaps
2571!> \param current_overlap current overlap in dbcsr format
2572!> \param coeffs resulting nvec coefficients
2573!> \param nvec number of previous overlaps
2574!> \param eps Tikhonov regularization
2575!> \param io_unit output unit
2576!> \param print_level print level
2577!> \param current_overlap_kp ...
2578!> \param kpoint_weights ...
2579!> \param para_env_inter_kp ...
2580!> \par History
2581!> 04.2026 created [Michele Nottoli]
2582! **************************************************************************************************
2583 SUBROUTINE tr_fitting(wf_history, current_overlap, coeffs, nvec, eps, io_unit, print_level, &
2584 current_overlap_kp, kpoint_weights, para_env_inter_kp)
2585 TYPE(qs_wf_history_type), POINTER :: wf_history
2586 TYPE(dbcsr_type), INTENT(IN) :: current_overlap
2587 INTEGER, INTENT(IN) :: nvec
2588 REAL(kind=dp), INTENT(OUT) :: coeffs(nvec)
2589 REAL(kind=dp), INTENT(IN) :: eps
2590 INTEGER, INTENT(IN) :: io_unit, print_level
2591 TYPE(cp_cfm_type), DIMENSION(:), INTENT(IN), &
2592 OPTIONAL :: current_overlap_kp
2593 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: kpoint_weights
2594 TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_inter_kp
2595
2596 COMPLEX(KIND=dp) :: ztrace
2597 INTEGER :: i, icol_local, ikp, info, irow_local, j, &
2598 ntr
2599 REAL(kind=dp) :: error, norm_ref, weight
2600 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: b
2601 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: a
2602 TYPE(cp_cfm_type) :: target_overlap_cfm, tmp_conj_cfm, &
2603 tmp_i_cfm, tmp_j_cfm
2604 TYPE(dbcsr_type) :: target_overlap, tmp_i, tmp_j, tmp_k
2605 TYPE(qs_wf_snapshot_type), POINTER :: ref_state, state
2606
2607 IF (nvec <= 0) THEN
2608 cpabort("Not enough vectors to do the fitting")
2609 ELSE IF (nvec == 1) THEN
2610 coeffs(1) = 1.0_dp
2611 RETURN
2612 END IF
2613
2614 IF (mod(nvec, 2) == 0) THEN
2615 ntr = nvec/2
2616 ELSE
2617 ntr = (nvec - 1)/2
2618 END IF
2619
2620 IF (PRESENT(current_overlap_kp)) THEN
2621 ALLOCATE (a(ntr, ntr), b(ntr))
2622 a = 0.0_dp
2623 b = 0.0_dp
2624
2625 ref_state => wfi_get_snapshot(wf_history, wf_index=nvec)
2626 CALL cp_cfm_create(target_overlap_cfm, current_overlap_kp(1)%matrix_struct)
2627 CALL cp_cfm_create(tmp_i_cfm, current_overlap_kp(1)%matrix_struct)
2628 CALL cp_cfm_create(tmp_j_cfm, current_overlap_kp(1)%matrix_struct)
2629 CALL cp_cfm_create(tmp_conj_cfm, current_overlap_kp(1)%matrix_struct)
2630
2631 DO ikp = 1, SIZE(current_overlap_kp)
2632 weight = 1.0_dp
2633 IF (PRESENT(kpoint_weights)) weight = kpoint_weights(ikp)
2634
2635 CALL cp_cfm_to_cfm(current_overlap_kp(ikp), target_overlap_cfm)
2636 CALL cp_cfm_scale_and_add(z_one, target_overlap_cfm, z_one, ref_state%overlap_cfm_kp(ikp))
2637 DO i = 1, ntr
2638 state => wfi_get_snapshot(wf_history, wf_index=i)
2639 CALL cp_cfm_to_cfm(state%overlap_cfm_kp(ikp), tmp_i_cfm)
2640 state => wfi_get_snapshot(wf_history, wf_index=nvec - i)
2641 CALL cp_cfm_scale_and_add(z_one, tmp_i_cfm, z_one, state%overlap_cfm_kp(ikp))
2642
2643 CALL cp_cfm_to_cfm(tmp_i_cfm, tmp_conj_cfm)
2644 DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
2645 DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
2646 tmp_conj_cfm%local_data(irow_local, icol_local) = &
2647 conjg(tmp_conj_cfm%local_data(irow_local, icol_local))
2648 END DO
2649 END DO
2650 CALL cp_cfm_trace(tmp_conj_cfm, target_overlap_cfm, ztrace)
2651 b(i) = b(i) + weight*real(ztrace, kind=dp)
2652 DO j = 1, i
2653 state => wfi_get_snapshot(wf_history, wf_index=j)
2654 CALL cp_cfm_to_cfm(state%overlap_cfm_kp(ikp), tmp_j_cfm)
2655 state => wfi_get_snapshot(wf_history, wf_index=nvec - j)
2656 CALL cp_cfm_scale_and_add(z_one, tmp_j_cfm, z_one, state%overlap_cfm_kp(ikp))
2657 CALL cp_cfm_trace(tmp_conj_cfm, tmp_j_cfm, ztrace)
2658 a(j, i) = a(j, i) + weight*real(ztrace, kind=dp)
2659 END DO
2660 END DO
2661 END DO
2662
2663 DO i = 1, ntr
2664 DO j = 1, i
2665 a(i, j) = a(j, i)
2666 END DO
2667 END DO
2668
2669 IF (PRESENT(para_env_inter_kp)) THEN
2670 IF (ASSOCIATED(para_env_inter_kp)) THEN
2671 CALL para_env_inter_kp%sum(a)
2672 CALL para_env_inter_kp%sum(b)
2673 END IF
2674 END IF
2675
2676 DO i = 1, ntr
2677 a(i, i) = a(i, i) + eps**2
2678 END DO
2679
2680 CALL dposv('u', ntr, 1, a, ntr, b, ntr, info)
2681 IF (info /= 0) THEN
2682 cpabort("DPOSV failed.")
2683 END IF
2684
2685 coeffs = 0.0_dp
2686 coeffs(nvec) = -1.0_dp
2687 DO i = 1, ntr
2688 coeffs(i) = coeffs(i) + b(i)
2689 coeffs(nvec - i) = coeffs(nvec - i) + b(i)
2690 END DO
2691
2692 IF (print_level > low_print_level) THEN
2693 error = 0.0_dp
2694 norm_ref = 0.0_dp
2695 DO ikp = 1, SIZE(current_overlap_kp)
2696 weight = 1.0_dp
2697 IF (PRESENT(kpoint_weights)) weight = kpoint_weights(ikp)
2698 CALL cp_cfm_to_cfm(current_overlap_kp(ikp), tmp_i_cfm)
2699 DO i = 1, nvec
2700 state => wfi_get_snapshot(wf_history, wf_index=i)
2701 CALL cp_cfm_scale_and_add(z_one, tmp_i_cfm, cmplx(-coeffs(i), 0.0_dp, kind=dp), &
2702 state%overlap_cfm_kp(ikp))
2703 END DO
2704 CALL cp_cfm_to_cfm(tmp_i_cfm, tmp_conj_cfm)
2705 DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
2706 DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
2707 tmp_conj_cfm%local_data(irow_local, icol_local) = &
2708 conjg(tmp_conj_cfm%local_data(irow_local, icol_local))
2709 END DO
2710 END DO
2711 CALL cp_cfm_trace(tmp_conj_cfm, tmp_i_cfm, ztrace)
2712 error = error + weight*real(ztrace, kind=dp)
2713 CALL cp_cfm_to_cfm(ref_state%overlap_cfm_kp(ikp), tmp_conj_cfm)
2714 DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
2715 DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
2716 tmp_conj_cfm%local_data(irow_local, icol_local) = &
2717 conjg(tmp_conj_cfm%local_data(irow_local, icol_local))
2718 END DO
2719 END DO
2720 CALL cp_cfm_trace(tmp_conj_cfm, ref_state%overlap_cfm_kp(ikp), ztrace)
2721 norm_ref = norm_ref + weight*real(ztrace, kind=dp)
2722 END DO
2723 IF (PRESENT(para_env_inter_kp)) THEN
2724 IF (ASSOCIATED(para_env_inter_kp)) THEN
2725 CALL para_env_inter_kp%sum(error)
2726 CALL para_env_inter_kp%sum(norm_ref)
2727 END IF
2728 END IF
2729 IF (io_unit > 0) THEN
2730 WRITE (unit=io_unit, fmt="(/,T2,A,F20.10)") "GEXT overlap fitting error:", &
2731 sqrt(error/max(norm_ref, tiny(1.0_dp)))
2732 END IF
2733 END IF
2734
2735 CALL cp_cfm_release(target_overlap_cfm)
2736 CALL cp_cfm_release(tmp_i_cfm)
2737 CALL cp_cfm_release(tmp_j_cfm)
2738 CALL cp_cfm_release(tmp_conj_cfm)
2739 DEALLOCATE (a, b)
2740 RETURN
2741 END IF
2742
2743 ALLOCATE (a(ntr, ntr), b(ntr))
2744
2745 ! get the reference for the difference fitting
2746 ref_state => wfi_get_snapshot(wf_history, wf_index=nvec)
2747
2748 ! assemble the target sum
2749 CALL dbcsr_copy(target_overlap, current_overlap)
2750 CALL dbcsr_add(target_overlap, ref_state%overlap, 1.0_dp, 1.0_dp)
2751
2752 ! allocate tmp_k
2753 CALL dbcsr_copy(tmp_k, current_overlap)
2754
2755 ! assemble the matrix A and the RHS b
2756 DO i = 1, ntr
2757 state => wfi_get_snapshot(wf_history, wf_index=i)
2758 CALL dbcsr_copy(tmp_i, state%overlap)
2759 state => wfi_get_snapshot(wf_history, wf_index=nvec - i)
2760 CALL dbcsr_add(tmp_i, state%overlap, 1.0_dp, 1.0_dp)
2761
2762 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_i, target_overlap, &
2763 0.0_dp, tmp_k)
2764 CALL dbcsr_trace(tmp_k, b(i))
2765 DO j = 1, i
2766 state => wfi_get_snapshot(wf_history, wf_index=j)
2767 CALL dbcsr_copy(tmp_j, state%overlap)
2768 state => wfi_get_snapshot(wf_history, wf_index=nvec - j)
2769 CALL dbcsr_add(tmp_j, state%overlap, 1.0_dp, 1.0_dp)
2770 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_i, tmp_j, 0.0_dp, tmp_k)
2771 CALL dbcsr_trace(tmp_k, a(j, i))
2772 a(i, j) = a(j, i)
2773 END DO
2774 END DO
2775
2776 ! add the Tikhonov regularization
2777 DO i = 1, ntr
2778 a(i, i) = a(i, i) + eps**2
2779 END DO
2780
2781 ! solve the linear system
2782 CALL dposv('u', ntr, 1, a, ntr, b, ntr, info)
2783 IF (info /= 0) THEN
2784 cpabort("DPOSV failed.")
2785 END IF
2786
2787 ! reorder the coefficients
2788 coeffs = 0.0_dp
2789 coeffs(nvec) = -1.0_dp
2790 DO i = 1, ntr
2791 coeffs(i) = coeffs(i) + b(i)
2792 coeffs(nvec - i) = coeffs(nvec - i) + b(i)
2793 END DO
2794
2795 ! as a consistency check, print how well the current overlap
2796 ! is approximated by the linear combination of previous overlaps
2797 IF (print_level > low_print_level) THEN
2798 CALL dbcsr_copy(tmp_i, current_overlap)
2799 DO i = 1, nvec
2800 state => wfi_get_snapshot(wf_history, wf_index=i)
2801 CALL dbcsr_add(tmp_i, state%overlap, 1.0_dp, -coeffs(i))
2802 END DO
2803 error = dbcsr_frobenius_norm(tmp_i)/dbcsr_frobenius_norm(state%overlap)
2804 IF (io_unit > 0) THEN
2805 WRITE (unit=io_unit, fmt="(/,T2,A,F20.10)") "GEXT overlap fitting error:", error
2806 END IF
2807 END IF
2808
2809 ! free the memory
2810 CALL dbcsr_release(tmp_i)
2811 CALL dbcsr_release(tmp_j)
2812 CALL dbcsr_release(tmp_k)
2813 CALL dbcsr_release(target_overlap)
2814 DEALLOCATE (a, b)
2815
2816 END SUBROUTINE tr_fitting
2817
2818END 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:89
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:84
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
pw_methods::pw_copy
Definition pw_methods.F:146
qs_density_matrices::calculate_density_matrix
Definition qs_density_matrices.F:58
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::vandevondele2005a
integer, save, public vandevondele2005a
Definition bibliography.F:36
bibliography::kuhne2007
integer, save, public kuhne2007
Definition bibliography.F:36
bibliography::kolafa2004
integer, save, public kolafa2004
Definition bibliography.F:36
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:563
cp_cfm_basic_linalg::cp_cfm_scale_and_add_fm
subroutine, public cp_cfm_scale_and_add_fm(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b). where b is a real matrix (adapted from cp_cf...
Definition cp_cfm_basic_linalg.F:351
cp_cfm_basic_linalg::cp_cfm_triangular_multiply
subroutine, public cp_cfm_triangular_multiply(triangular_matrix, matrix_b, side, transa_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, alpha)
Multiplies in place by a triangular matrix: matrix_b = alpha op(triangular_matrix) matrix_b or (if si...
Definition cp_cfm_basic_linalg.F:980
cp_cfm_basic_linalg::cp_cfm_column_scale
subroutine, public cp_cfm_column_scale(matrix_a, scaling)
Scales columns of the full matrix by corresponding factors.
Definition cp_cfm_basic_linalg.F:633
cp_cfm_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:910
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:840
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:125
cp_cfm_types::cp_cfm_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:630
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:493
cp_cfm_types::cp_cfm_set_all
subroutine, public cp_cfm_set_all(matrix, alpha, beta)
Set all elements of the full matrix to alpha. Besides, set all diagonal matrix elements to beta (if g...
Definition cp_cfm_types.F:202
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition cp_cfm_types.F:788
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:524
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:370
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1095
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:807
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1203
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1141
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:253
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_trace
subroutine, public dbcsr_trace(matrix, trace)
Computes the trace of the given matrix, also known as the sum of its diagonal elements.
Definition cp_dbcsr_contrib.F:412
cp_dbcsr_contrib::dbcsr_frobenius_norm
real(dp) function, public dbcsr_frobenius_norm(matrix)
Compute the frobenius norm of a dbcsr matrix.
Definition cp_dbcsr_contrib.F:127
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:552
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:214
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:172
cp_fm_basic_linalg::cp_fm_scale
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
Definition cp_fm_basic_linalg.F:1730
cp_fm_pool_types
pool for for elements that are retained and released
Definition cp_fm_pool_types.F:14
cp_fm_pool_types::fm_pool_create_fm
subroutine, public fm_pool_create_fm(pool, element, name)
returns an element, allocating it if none is in the pool
Definition cp_fm_pool_types.F:185
cp_fm_pool_types::fm_pool_give_back_fm
subroutine, public fm_pool_give_back_fm(pool, element)
returns the element to the pool
Definition cp_fm_pool_types.F:224
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_start_copy_general
subroutine, public cp_fm_start_copy_general(source, destination, para_env, info)
Initiates the copy operation: get distribution data, post MPI isend and irecvs.
Definition cp_fm_types.F:1645
cp_fm_types::cp_fm_cleanup_copy_general
subroutine, public cp_fm_cleanup_copy_general(info)
Completes the copy operation: wait for comms clean up MPI state.
Definition cp_fm_types.F:2023
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1087
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:169
cp_fm_types::cp_fm_set_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:789
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:556
cp_fm_types::cp_fm_finish_copy_general
subroutine, public cp_fm_finish_copy_general(destination, info)
Completes the copy operation: wait for comms, unpack, clean up MPI state.
Definition cp_fm_types.F:1959
cp_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:972
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:1067
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::wfi_frozen_method_nr
integer, parameter, public wfi_frozen_method_nr
Definition input_constants.F:246
input_constants::wfi_linear_wf_method_nr
integer, parameter, public wfi_linear_wf_method_nr
Definition input_constants.F:246
input_constants::wfi_linear_p_method_nr
integer, parameter, public wfi_linear_p_method_nr
Definition input_constants.F:246
input_constants::wfi_linear_ps_method_nr
integer, parameter, public wfi_linear_ps_method_nr
Definition input_constants.F:246
input_constants::wfi_use_prev_rho_r_method_nr
integer, parameter, public wfi_use_prev_rho_r_method_nr
Definition input_constants.F:246
input_constants::wfi_use_guess_method_nr
integer, parameter, public wfi_use_guess_method_nr
Definition input_constants.F:246
input_constants::wfi_gext_proj_qtr_nr
integer, parameter, public wfi_gext_proj_qtr_nr
Definition input_constants.F:246
input_constants::wfi_use_prev_wf_method_nr
integer, parameter, public wfi_use_prev_wf_method_nr
Definition input_constants.F:246
input_constants::wfi_ps_method_nr
integer, parameter, public wfi_ps_method_nr
Definition input_constants.F:246
input_constants::wfi_gext_proj_nr
integer, parameter, public wfi_gext_proj_nr
Definition input_constants.F:246
input_constants::wfi_aspc_nr
integer, parameter, public wfi_aspc_nr
Definition input_constants.F:246
input_constants::wfi_use_prev_p_method_nr
integer, parameter, public wfi_use_prev_p_method_nr
Definition input_constants.F:246
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kpoint_methods
Routines needed for kpoint calculation.
Definition kpoint_methods.F:15
kpoint_methods::rskp_transform
subroutine, public rskp_transform(rmatrix, cmatrix, rsmat, ispin, xkp, cell_to_index, sab_nl, is_complex, rs_sign)
Transformation of real space matrices to a kpoint.
Definition kpoint_methods.F:1114
kpoint_methods::kpoint_density_transform
subroutine, public kpoint_density_transform(kpoint, denmat, wtype, tempmat, sab_nl, fmwork, for_aux_fit, pmat_ext)
generate real space density matrices in DBCSR format
Definition kpoint_methods.F:1741
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:1465
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:1214
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:525
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:213
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, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:530
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, mimic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, rhoz_cneo_set, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, harris_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Set the QUICKSTEP environment.
Definition qs_environment_types.F:1106
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::qs_ks_did_change
subroutine, public qs_ks_did_change(ks_env, s_mstruct_changed, rho_changed, potential_changed, full_reset)
tells that some of the things relevant to the ks calculation did change. has to be called when change...
Definition qs_ks_types.F:956
qs_matrix_pools
wrapper for the pools of matrixes
Definition qs_matrix_pools.F:14
qs_matrix_pools::mpools_get
subroutine, public mpools_get(mpools, ao_mo_fm_pools, ao_ao_fm_pools, mo_mo_fm_pools, ao_mosub_fm_pools, mosub_mosub_fm_pools, maxao_maxmo_fm_pool, maxao_maxao_fm_pool, maxmo_maxmo_fm_pool)
returns various attributes of the mpools (notably the pools contained in it)
Definition qs_matrix_pools.F:139
qs_mo_methods
collects routines that perform operations directly related to MOs
Definition qs_mo_methods.F:14
qs_mo_methods::make_basis_simple
subroutine, public make_basis_simple(vmatrix, ncol)
given a set of vectors, return an orthogonal (C^T C == 1) set spanning the same space (notice,...
Definition qs_mo_methods.F:356
qs_mo_methods::make_basis_lowdin
subroutine, public make_basis_lowdin(vmatrix, ncol, matrix_s)
return a set of S orthonormal vectors (C^T S C == 1) where a Loedwin transformation is applied to kee...
Definition qs_mo_methods.F:300
qs_mo_methods::make_basis_cholesky
subroutine, public make_basis_cholesky(vmatrix, ncol, ortho)
return a set of S orthonormal vectors (C^T S C == 1) where the cholesky decomposed form of S is passe...
Definition qs_mo_methods.F:249
qs_mo_methods::make_basis_sm
subroutine, public make_basis_sm(vmatrix, ncol, matrix_s)
returns an S-orthonormal basis v (v^T S v ==1)
Definition qs_mo_methods.F:89
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F: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:378
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_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:2277
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:2210
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:2182
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:744
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_pool_types::cp_fm_pool_p_type
to create arrays of pools
Definition cp_fm_pool_types.F:70
cp_fm_pool_types::cp_fm_pool_type
represent a pool of elements with the same structure
Definition cp_fm_pool_types.F:56
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::copy_info_type
Stores the state of a copy between cp_fm_start_copy_general and cp_fm_finish_copy_general.
Definition cp_fm_types.F:141
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:115
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
kpoint_types::kpoint_env_type
Keeps information about a specific k-point.
Definition kpoint_types.F:83
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:171
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
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:220
qs_matrix_pools::qs_matrix_pools_type
container for the pools of matrixes used by qs
Definition qs_matrix_pools.F:66
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:97
qs_wf_history_types::qs_wf_history_type
keeps track of the previous wavefunctions and can extrapolate them for the next step of md
Definition qs_wf_history_types.F: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