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