Line data Source code
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 : ! **************************************************************************************************
21 : MODULE 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 :
134 : CONTAINS
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 15034 : 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 15034 : snapshot%dt = 1.0_dp
152 15034 : 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 25480 : 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 25480 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
178 25480 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
179 : LOGICAL :: my_kpgrp
180 25480 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
181 : TYPE(cell_type), POINTER :: cell
182 25480 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info_ft
183 25480 : 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 25480 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
188 25480 : 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 25480 : 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 25480 : POINTER :: sab_nl
197 25480 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
198 25480 : 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 25480 : 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 25480 : CALL timeset(routineN, handle)
207 :
208 25480 : NULLIFY (pw_env, auxbas_pw_pool, ao_mo_pools, ao_ao_fm_pools, dft_control, mos, mo_coeff, &
209 25480 : rho, rho_r, rho_g, rho_ao, matrix_s, matrix_s_kp, kpoints, kp, cell, &
210 25480 : particle_set, kp_dist, cell_to_index, xkp, sab_nl, scf_env, mpools_kp, para_env_ft, &
211 25480 : rmat_ft, cmat_ft, tmpmat_ft, ao_ao_struct_ft)
212 : CALL get_qs_env(qs_env, pw_env=pw_env, &
213 25480 : dft_control=dft_control, rho=rho, cell=cell, particle_set=particle_set)
214 25480 : CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_pools)
215 25480 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
216 :
217 25480 : CPASSERT(ASSOCIATED(wf_history))
218 25480 : CPASSERT(ASSOCIATED(dft_control))
219 25480 : IF (.NOT. ASSOCIATED(snapshot)) THEN
220 15034 : ALLOCATE (snapshot)
221 15034 : CALL wfs_create(snapshot)
222 : END IF
223 25480 : CPASSERT(wf_history%ref_count > 0)
224 :
225 25480 : nspins = dft_control%nspins
226 25480 : snapshot%dt = 1.0_dp
227 25480 : IF (PRESENT(dt)) snapshot%dt = dt
228 25480 : IF (wf_history%store_wf) THEN
229 21980 : CALL get_qs_env(qs_env, mos=mos)
230 21980 : IF (.NOT. ASSOCIATED(snapshot%wf)) THEN
231 : CALL fm_pools_create_fm_vect(ao_mo_pools, snapshot%wf, &
232 12328 : name="ws_snap-ws")
233 12328 : CPASSERT(nspins == SIZE(snapshot%wf))
234 : END IF
235 46976 : DO ispin = 1, nspins
236 24996 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
237 46976 : CALL cp_fm_to_fm(mo_coeff, snapshot%wf(ispin))
238 : END DO
239 : ELSE
240 3500 : CALL fm_pools_give_back_fm_vect(ao_mo_pools, snapshot%wf)
241 : END IF
242 :
243 25480 : IF (wf_history%store_rho_r) THEN
244 0 : CALL qs_rho_get(rho, rho_r=rho_r)
245 0 : CPASSERT(ASSOCIATED(rho_r))
246 0 : IF (.NOT. ASSOCIATED(snapshot%rho_r)) THEN
247 0 : ALLOCATE (snapshot%rho_r(nspins))
248 0 : DO ispin = 1, nspins
249 0 : CALL auxbas_pw_pool%create_pw(snapshot%rho_r(ispin))
250 : END DO
251 : END IF
252 0 : DO ispin = 1, nspins
253 0 : CALL pw_copy(rho_r(ispin), snapshot%rho_r(ispin))
254 : END DO
255 25480 : ELSE IF (ASSOCIATED(snapshot%rho_r)) THEN
256 0 : DO ispin = 1, SIZE(snapshot%rho_r)
257 0 : CALL auxbas_pw_pool%give_back_pw(snapshot%rho_r(ispin))
258 : END DO
259 0 : DEALLOCATE (snapshot%rho_r)
260 : END IF
261 :
262 25480 : IF (wf_history%store_rho_g) THEN
263 0 : CALL qs_rho_get(rho, rho_g=rho_g)
264 0 : CPASSERT(ASSOCIATED(rho_g))
265 0 : IF (.NOT. ASSOCIATED(snapshot%rho_g)) THEN
266 0 : ALLOCATE (snapshot%rho_g(nspins))
267 0 : DO ispin = 1, nspins
268 0 : CALL auxbas_pw_pool%create_pw(snapshot%rho_g(ispin))
269 : END DO
270 : END IF
271 0 : DO ispin = 1, nspins
272 0 : CALL pw_copy(rho_g(ispin), snapshot%rho_g(ispin))
273 : END DO
274 25480 : ELSE IF (ASSOCIATED(snapshot%rho_g)) THEN
275 0 : DO ispin = 1, SIZE(snapshot%rho_g)
276 0 : CALL auxbas_pw_pool%give_back_pw(snapshot%rho_g(ispin))
277 : END DO
278 0 : DEALLOCATE (snapshot%rho_g)
279 : END IF
280 :
281 25480 : IF (ASSOCIATED(snapshot%rho_ao)) THEN ! the sparsity might be different
282 : ! (future struct:check)
283 262 : CALL dbcsr_deallocate_matrix_set(snapshot%rho_ao)
284 : END IF
285 25480 : IF (wf_history%store_rho_ao) THEN
286 326 : CALL qs_rho_get(rho, rho_ao=rho_ao)
287 326 : CPASSERT(ASSOCIATED(rho_ao))
288 :
289 326 : CALL dbcsr_allocate_matrix_set(snapshot%rho_ao, nspins)
290 812 : DO ispin = 1, nspins
291 486 : ALLOCATE (snapshot%rho_ao(ispin)%matrix)
292 812 : CALL dbcsr_copy(snapshot%rho_ao(ispin)%matrix, rho_ao(ispin)%matrix)
293 : END DO
294 : END IF
295 :
296 25480 : IF (ASSOCIATED(snapshot%rho_ao_kp)) THEN ! the sparsity might be different
297 : ! (future struct:check)
298 220 : CALL dbcsr_deallocate_matrix_set(snapshot%rho_ao_kp)
299 : END IF
300 25480 : IF (wf_history%store_rho_ao_kp) THEN
301 232 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
302 232 : CPASSERT(ASSOCIATED(rho_ao_kp))
303 :
304 232 : nimg = dft_control%nimages
305 232 : CALL dbcsr_allocate_matrix_set(snapshot%rho_ao_kp, nspins, nimg)
306 554 : DO ispin = 1, nspins
307 34092 : DO img = 1, nimg
308 33538 : ALLOCATE (snapshot%rho_ao_kp(ispin, img)%matrix)
309 : CALL dbcsr_copy(snapshot%rho_ao_kp(ispin, img)%matrix, &
310 33860 : rho_ao_kp(ispin, img)%matrix)
311 : END DO
312 : END DO
313 : END IF
314 :
315 25480 : IF (ASSOCIATED(snapshot%overlap)) THEN ! the sparsity might be different
316 : ! (future struct:check)
317 7044 : CALL dbcsr_deallocate_matrix(snapshot%overlap)
318 : END IF
319 25480 : IF (wf_history%store_overlap) THEN
320 18338 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
321 18338 : CPASSERT(ASSOCIATED(matrix_s))
322 18338 : CPASSERT(ASSOCIATED(matrix_s(1)%matrix))
323 18338 : ALLOCATE (snapshot%overlap)
324 18338 : CALL dbcsr_copy(snapshot%overlap, matrix_s(1)%matrix)
325 : END IF
326 :
327 25480 : CALL get_qs_env(qs_env, kpoints=kpoints)
328 25480 : IF (ASSOCIATED(kpoints)) THEN
329 25480 : IF (ASSOCIATED(kpoints%kp_env)) THEN
330 : ! --- k-point WFN snapshot: store complex MO coefficients per local k-point ---
331 3158 : IF (wf_history%store_wf_kp) THEN
332 2926 : CALL get_kpoint_info(kpoints, kp_range=kp_range)
333 2926 : kplocal = kp_range(2) - kp_range(1) + 1
334 2926 : nspin_kp = SIZE(kpoints%kp_env(1)%kpoint_env%mos, 2)
335 2926 : nc = SIZE(kpoints%kp_env(1)%kpoint_env%mos, 1) ! 2=complex, 1=real
336 :
337 2926 : CALL wfi_store_kp_pbc_shift(snapshot, cell, particle_set)
338 :
339 2926 : IF (ASSOCIATED(snapshot%wf_kp)) THEN
340 842 : DO ikp = 1, SIZE(snapshot%wf_kp, 1)
341 1926 : DO ic = 1, SIZE(snapshot%wf_kp, 2)
342 2710 : DO ispin = 1, SIZE(snapshot%wf_kp, 3)
343 2168 : CALL cp_fm_release(snapshot%wf_kp(ikp, ic, ispin))
344 : END DO
345 : END DO
346 : END DO
347 300 : DEALLOCATE (snapshot%wf_kp)
348 : END IF
349 :
350 36356 : ALLOCATE (snapshot%wf_kp(kplocal, nc, nspin_kp))
351 8410 : DO ikp = 1, kplocal
352 5484 : kp => kpoints%kp_env(ikp)%kpoint_env
353 14628 : DO ispin = 1, nspin_kp
354 24136 : DO ic = 1, nc
355 12434 : 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 12434 : name="wfkp_snap")
359 18652 : 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 3158 : IF (wf_history%store_overlap_kp) THEN
370 2888 : 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 2888 : sab_nl=sab_nl, cell_to_index=cell_to_index)
374 2888 : kplocal = kp_range(2) - kp_range(1) + 1
375 2888 : 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 2888 : ALLOCATE (rmat_ft, cmat_ft, tmpmat_ft)
379 : CALL dbcsr_create(rmat_ft, template=matrix_s_kp(1, 1)%matrix, &
380 2888 : matrix_type=dbcsr_type_symmetric)
381 : CALL dbcsr_create(cmat_ft, template=matrix_s_kp(1, 1)%matrix, &
382 2888 : matrix_type=dbcsr_type_antisymmetric)
383 : CALL dbcsr_create(tmpmat_ft, template=matrix_s_kp(1, 1)%matrix, &
384 2888 : matrix_type=dbcsr_type_no_symmetry)
385 2888 : CALL cp_dbcsr_alloc_block_from_nbl(rmat_ft, sab_nl)
386 2888 : CALL cp_dbcsr_alloc_block_from_nbl(cmat_ft, sab_nl)
387 :
388 : ! Get kp-subgroup FM from pool
389 2888 : CALL get_kpoint_info(kpoints, mpools=mpools_kp)
390 2888 : CALL mpools_get(mpools_kp, ao_ao_fm_pools=ao_ao_fm_pools)
391 2888 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal_ft)
392 :
393 : ! Release old snapshot if present
394 2888 : IF (ASSOCIATED(snapshot%overlap_cfm_kp)) THEN
395 822 : DO ikp = 1, SIZE(snapshot%overlap_cfm_kp)
396 822 : CALL cp_cfm_release(snapshot%overlap_cfm_kp(ikp))
397 : END DO
398 296 : DEALLOCATE (snapshot%overlap_cfm_kp)
399 : END IF
400 14084 : ALLOCATE (snapshot%overlap_cfm_kp(kplocal))
401 :
402 2888 : CALL cp_fm_get_info(fmlocal_ft, matrix_struct=ao_ao_struct_ft)
403 :
404 : ! Communication info array
405 55884 : ALLOCATE (info_ft(kplocal*nkp_grps, 2))
406 :
407 : ! Phase A: Start async FT + redistribute for each k-point
408 2888 : indx_ft = 0
409 8308 : DO ikp = 1, kplocal
410 17478 : DO igroup = 1, nkp_grps
411 9170 : ik = kp_dist(1, igroup) + ikp - 1
412 9170 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
413 9170 : indx_ft = indx_ft + 1
414 :
415 9170 : CALL dbcsr_set(rmat_ft, 0.0_dp)
416 9170 : 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 9170 : cell_to_index=cell_to_index, sab_nl=sab_nl)
420 9170 : CALL dbcsr_desymmetrize(rmat_ft, tmpmat_ft)
421 9170 : CALL copy_dbcsr_to_fm(tmpmat_ft, scf_env%scf_work1(1))
422 9170 : CALL dbcsr_desymmetrize(cmat_ft, tmpmat_ft)
423 9170 : CALL copy_dbcsr_to_fm(tmpmat_ft, scf_env%scf_work1(2))
424 :
425 14590 : IF (my_kpgrp) THEN
426 : CALL cp_fm_start_copy_general(scf_env%scf_work1(1), fmlocal_ft, &
427 5420 : para_env_ft, info_ft(indx_ft, 1))
428 : CALL cp_fm_start_copy_general(scf_env%scf_work1(2), fmlocal_ft, &
429 5420 : 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 3750 : para_env_ft, info_ft(indx_ft, 1))
433 : CALL cp_fm_start_copy_general(scf_env%scf_work1(2), fmdummy_ft, &
434 3750 : 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 8308 : DO ikp = 1, kplocal
442 5420 : CALL cp_cfm_create(snapshot%overlap_cfm_kp(ikp), ao_ao_struct_ft)
443 5420 : CALL cp_cfm_set_all(snapshot%overlap_cfm_kp(ikp), z_zero)
444 17478 : DO igroup = 1, nkp_grps
445 9170 : ik = kp_dist(1, igroup) + ikp - 1
446 9170 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
447 3750 : indx_ft = indx_ft + 1
448 5420 : IF (my_kpgrp) THEN
449 5420 : 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 5420 : z_one, fmlocal_ft)
452 5420 : 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 5420 : gaussi, fmlocal_ft)
455 : END IF
456 : END DO
457 : END DO
458 :
459 : ! Cleanup
460 12058 : DO indx_ft = 1, kplocal*nkp_grps
461 9170 : CALL cp_fm_cleanup_copy_general(info_ft(indx_ft, 1))
462 12058 : CALL cp_fm_cleanup_copy_general(info_ft(indx_ft, 2))
463 : END DO
464 24116 : DEALLOCATE (info_ft)
465 2888 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal_ft)
466 2888 : CALL dbcsr_deallocate_matrix(rmat_ft)
467 2888 : CALL dbcsr_deallocate_matrix(cmat_ft)
468 5776 : 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 25480 : CALL timestop(handle)
479 :
480 25480 : 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 9066 : 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 9066 : CALL timeset(routineN, handle)
506 :
507 9066 : ALLOCATE (wf_history)
508 9066 : wf_history%ref_count = 1
509 9066 : wf_history%memory_depth = 0
510 9066 : wf_history%snapshot_count = 0
511 9066 : 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 9066 : 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 64 : wf_history%memory_depth = 0
530 : CASE (wfi_use_prev_p_method_nr)
531 64 : wf_history%memory_depth = 1
532 64 : wf_history%store_rho_ao = .TRUE.
533 : CASE (wfi_linear_wf_method_nr)
534 4 : wf_history%memory_depth = 2
535 4 : wf_history%store_wf = .TRUE.
536 : CASE (wfi_linear_p_method_nr)
537 6 : wf_history%memory_depth = 2
538 6 : wf_history%store_rho_ao = .TRUE.
539 : CASE (wfi_linear_ps_method_nr)
540 6 : wf_history%memory_depth = 2
541 6 : wf_history%store_wf = .TRUE.
542 6 : IF (.NOT. has_unit_metric) wf_history%store_overlap = .TRUE.
543 : CASE (wfi_ps_method_nr)
544 345 : CALL cite_reference(VandeVondele2005a)
545 345 : wf_history%memory_depth = extrapolation_order + 1
546 345 : wf_history%store_wf = .TRUE.
547 345 : wf_history%store_wf_kp = .TRUE.
548 345 : IF (.NOT. has_unit_metric) THEN
549 341 : wf_history%store_overlap = .TRUE.
550 341 : wf_history%store_overlap_kp = .TRUE.
551 : END IF
552 : CASE (wfi_frozen_method_nr)
553 4 : wf_history%memory_depth = 1
554 4 : wf_history%store_frozen_density = .TRUE.
555 : CASE (wfi_aspc_nr)
556 8265 : wf_history%memory_depth = extrapolation_order + 2
557 8265 : wf_history%store_wf = .TRUE.
558 8265 : wf_history%store_wf_kp = .TRUE.
559 8265 : IF (.NOT. has_unit_metric) THEN
560 7283 : wf_history%store_overlap = .TRUE.
561 7283 : wf_history%store_overlap_kp = .TRUE.
562 : END IF
563 : CASE (wfi_gext_proj_nr)
564 24 : wf_history%memory_depth = extrapolation_order
565 24 : wf_history%store_wf = .TRUE.
566 24 : wf_history%store_wf_kp = .TRUE.
567 24 : wf_history%store_overlap = .TRUE.
568 24 : wf_history%store_overlap_kp = .TRUE.
569 : CASE (wfi_gext_proj_qtr_nr)
570 6 : wf_history%memory_depth = extrapolation_order
571 6 : wf_history%store_wf = .TRUE.
572 6 : wf_history%store_wf_kp = .TRUE.
573 6 : wf_history%store_overlap = .TRUE.
574 6 : wf_history%store_overlap_kp = .TRUE.
575 : CASE default
576 : CALL cp_abort(__LOCATION__, &
577 : "Unknown interpolation method: "// &
578 9066 : TRIM(ADJUSTL(cp_to_string(interpolation_method_nr))))
579 : END SELECT
580 69445 : ALLOCATE (wf_history%past_states(wf_history%memory_depth))
581 :
582 51655 : DO i = 1, SIZE(wf_history%past_states)
583 51655 : NULLIFY (wf_history%past_states(i)%snapshot)
584 : END DO
585 :
586 9066 : CALL timestop(handle)
587 9066 : 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 654 : SUBROUTINE wfi_create_for_kp(wf_history)
599 : TYPE(qs_wf_history_type), POINTER :: wf_history
600 :
601 : INTEGER :: i
602 :
603 654 : CPASSERT(ASSOCIATED(wf_history))
604 654 : IF (wf_history%store_rho_ao) THEN
605 10 : wf_history%store_rho_ao_kp = .TRUE.
606 10 : 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 654 : IF (wf_history%interpolation_method_nr == wfi_use_prev_wf_method_nr) THEN
612 30 : wf_history%memory_depth = 1
613 30 : wf_history%store_wf_kp = .TRUE.
614 30 : wf_history%store_wf = .FALSE.
615 30 : wf_history%store_overlap = .FALSE.
616 30 : IF (ASSOCIATED(wf_history%past_states)) DEALLOCATE (wf_history%past_states)
617 120 : ALLOCATE (wf_history%past_states(wf_history%memory_depth))
618 60 : DO i = 1, SIZE(wf_history%past_states)
619 60 : NULLIFY (wf_history%past_states(i)%snapshot)
620 : END DO
621 624 : ELSE IF (wf_history%store_wf_kp) THEN
622 478 : wf_history%store_wf = .FALSE.
623 478 : 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 146 : IF (wf_history%store_wf .OR. wf_history%store_overlap) THEN
628 0 : CPABORT("Linear WFN-based extrapolation methods not implemented for k-points.")
629 : END IF
630 : END IF
631 654 : IF (wf_history%store_frozen_density) THEN
632 0 : CPABORT("Frozen density initialization method not possible for kpoints.")
633 : END IF
634 :
635 654 : 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 13485 : FUNCTION wfi_get_method_label(method_nr) RESULT(res)
646 : INTEGER, INTENT(in) :: method_nr
647 : CHARACTER(len=30) :: res
648 :
649 13485 : res = "unknown"
650 13723 : SELECT CASE (method_nr)
651 : CASE (wfi_use_prev_p_method_nr)
652 238 : res = "previous_p"
653 : CASE (wfi_use_prev_wf_method_nr)
654 237 : res = "previous_wf"
655 : CASE (wfi_use_guess_method_nr)
656 5138 : res = "initial_guess"
657 : CASE (wfi_linear_wf_method_nr)
658 2 : res = "mo linear"
659 : CASE (wfi_linear_p_method_nr)
660 3 : res = "P linear"
661 : CASE (wfi_linear_ps_method_nr)
662 6 : res = "PS linear"
663 : CASE (wfi_ps_method_nr)
664 188 : res = "PS Nth order"
665 : CASE (wfi_frozen_method_nr)
666 4 : res = "frozen density approximation"
667 : CASE (wfi_aspc_nr)
668 7584 : res = "ASPC"
669 : CASE (wfi_gext_proj_nr)
670 71 : res = "GEXT_PROJ"
671 : CASE (wfi_gext_proj_qtr_nr)
672 14 : res = "GEXT_PROJ_QTR"
673 : CASE default
674 : CALL cp_abort(__LOCATION__, &
675 : "Unknown interpolation method: "// &
676 13485 : TRIM(ADJUSTL(cp_to_string(method_nr))))
677 : END SELECT
678 13485 : 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 26671 : 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 26671 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coeffs
711 26671 : 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 26671 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao, rho_frozen_ao
717 26671 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
718 26671 : 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 26671 : NULLIFY (mos, ao_mo_fm_pools, t0_state, t1_state, mo_coeff, &
723 26671 : rho, rho_ao, rho_frozen_ao)
724 :
725 26671 : use_overlap = wf_history%store_overlap
726 :
727 26671 : CALL timeset(routineN, handle)
728 26671 : logger => cp_get_default_logger()
729 26671 : 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 26671 : extension=".scfLog")
732 :
733 26671 : CPASSERT(ASSOCIATED(wf_history))
734 26671 : CPASSERT(wf_history%ref_count > 0)
735 26671 : CPASSERT(ASSOCIATED(qs_env))
736 26671 : CALL get_qs_env(qs_env, mos=mos, rho=rho, do_kpoints=do_kpoints)
737 26671 : CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
738 : ! chooses the method for this extrapolation
739 26671 : IF (wf_history%snapshot_count < 1) THEN
740 : actual_extrapolation_method_nr = wfi_use_guess_method_nr
741 : ELSE
742 16560 : actual_extrapolation_method_nr = wf_history%interpolation_method_nr
743 : END IF
744 :
745 8 : SELECT CASE (actual_extrapolation_method_nr)
746 : CASE (wfi_linear_wf_method_nr)
747 8 : IF (wf_history%snapshot_count < 2) THEN
748 4 : actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
749 : END IF
750 : CASE (wfi_linear_p_method_nr)
751 12 : IF (wf_history%snapshot_count < 2) THEN
752 6 : actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
753 : END IF
754 : CASE (wfi_linear_ps_method_nr)
755 16560 : IF (wf_history%snapshot_count < 2) THEN
756 6 : actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
757 : END IF
758 : END SELECT
759 :
760 26671 : IF (PRESENT(extrapolation_method_nr)) THEN
761 26671 : extrapolation_method_nr = actual_extrapolation_method_nr
762 : END IF
763 26671 : my_orthogonal_wf = .FALSE.
764 :
765 8 : SELECT CASE (actual_extrapolation_method_nr)
766 : CASE (wfi_frozen_method_nr)
767 8 : CPASSERT(.NOT. do_kpoints)
768 8 : t0_state => wfi_get_snapshot(wf_history, wf_index=1)
769 8 : CPASSERT(ASSOCIATED(t0_state%rho_frozen))
770 :
771 8 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
772 8 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
773 :
774 8 : CALL qs_rho_get(t0_state%rho_frozen, rho_ao=rho_frozen_ao)
775 8 : CALL qs_rho_get(rho, rho_ao=rho_ao)
776 16 : DO ispin = 1, SIZE(rho_frozen_ao)
777 : CALL dbcsr_copy(rho_ao(ispin)%matrix, &
778 : rho_frozen_ao(ispin)%matrix, &
779 16 : 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 8 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
784 8 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
785 :
786 8 : my_orthogonal_wf = .FALSE.
787 : CASE (wfi_use_prev_p_method_nr)
788 476 : t0_state => wfi_get_snapshot(wf_history, wf_index=1)
789 476 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
790 476 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
791 476 : IF (do_kpoints) THEN
792 218 : CPASSERT(ASSOCIATED(t0_state%rho_ao_kp))
793 218 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
794 524 : DO ispin = 1, SIZE(t0_state%rho_ao_kp, 1)
795 31248 : DO img = 1, SIZE(t0_state%rho_ao_kp, 2)
796 31030 : IF (img > SIZE(rho_ao_kp, 2)) THEN
797 18 : 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 30706 : keep_sparsity=.TRUE.)
802 : END IF
803 : END DO
804 : END DO
805 : ELSE
806 258 : CPASSERT(ASSOCIATED(t0_state%rho_ao))
807 258 : CALL qs_rho_get(rho, rho_ao=rho_ao)
808 646 : DO ispin = 1, SIZE(t0_state%rho_ao)
809 : CALL dbcsr_copy(rho_ao(ispin)%matrix, &
810 : t0_state%rho_ao(ispin)%matrix, &
811 646 : 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 476 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
817 476 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
818 : CASE (wfi_use_prev_wf_method_nr)
819 474 : my_orthogonal_wf = .TRUE.
820 474 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
821 474 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
822 :
823 474 : IF (do_kpoints) THEN
824 6 : CALL wfi_use_prev_wf_kp(qs_env, io_unit, print_level)
825 : ELSE
826 468 : CALL qs_rho_get(rho, rho_ao=rho_ao)
827 988 : DO ispin = 1, SIZE(mos)
828 520 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
829 520 : CALL reorthogonalize_vectors(qs_env, v_matrix=mo_coeff, n_col=nmo)
830 1508 : CALL calculate_density_matrix(mo_set=mos(ispin), density_matrix=rho_ao(ispin)%matrix)
831 : END DO
832 468 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
833 468 : 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 10257 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
843 10257 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
844 : CASE (wfi_linear_wf_method_nr)
845 4 : CPASSERT(.NOT. do_kpoints)
846 4 : t0_state => wfi_get_snapshot(wf_history, wf_index=2)
847 4 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
848 4 : CPASSERT(ASSOCIATED(t0_state))
849 4 : CPASSERT(ASSOCIATED(t1_state))
850 4 : CPASSERT(ASSOCIATED(t0_state%wf))
851 4 : CPASSERT(ASSOCIATED(t1_state%wf))
852 4 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
853 4 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
854 :
855 4 : my_orthogonal_wf = .TRUE.
856 4 : t0 = 0.0_dp
857 4 : t1 = t1_state%dt
858 4 : t2 = t1 + dt
859 4 : CALL qs_rho_get(rho, rho_ao=rho_ao)
860 8 : DO ispin = 1, SIZE(mos)
861 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, &
862 4 : 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 4 : 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 4 : beta=(t1 - t2)/(t1 - t0), matrix_b=t0_state%wf(ispin))
871 : CALL reorthogonalize_vectors(qs_env, &
872 : v_matrix=mo_coeff, &
873 4 : n_col=nmo)
874 : CALL calculate_density_matrix(mo_set=mos(ispin), &
875 12 : density_matrix=rho_ao(ispin)%matrix)
876 : END DO
877 4 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
878 :
879 : CALL qs_ks_did_change(qs_env%ks_env, &
880 4 : rho_changed=.TRUE.)
881 : CASE (wfi_linear_p_method_nr)
882 6 : t0_state => wfi_get_snapshot(wf_history, wf_index=2)
883 6 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
884 6 : CPASSERT(ASSOCIATED(t0_state))
885 6 : CPASSERT(ASSOCIATED(t1_state))
886 6 : IF (do_kpoints) THEN
887 2 : CPASSERT(ASSOCIATED(t0_state%rho_ao_kp))
888 2 : CPASSERT(ASSOCIATED(t1_state%rho_ao_kp))
889 : ELSE
890 4 : CPASSERT(ASSOCIATED(t0_state%rho_ao))
891 4 : CPASSERT(ASSOCIATED(t1_state%rho_ao))
892 : END IF
893 6 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
894 6 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
895 :
896 6 : t0 = 0.0_dp
897 6 : t1 = t1_state%dt
898 6 : t2 = t1 + dt
899 6 : IF (do_kpoints) THEN
900 2 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
901 4 : DO ispin = 1, SIZE(rho_ao_kp, 1)
902 528 : DO img = 1, SIZE(rho_ao_kp, 2)
903 524 : IF (img > SIZE(t0_state%rho_ao_kp, 2) .OR. &
904 2 : img > SIZE(t1_state%rho_ao_kp, 2)) THEN
905 22 : 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 502 : 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 502 : alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0))
911 : END IF
912 : END DO
913 : END DO
914 : ELSE
915 4 : CALL qs_rho_get(rho, rho_ao=rho_ao)
916 8 : DO ispin = 1, SIZE(rho_ao)
917 : CALL dbcsr_add(rho_ao(ispin)%matrix, t1_state%rho_ao(ispin)%matrix, &
918 4 : 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 8 : alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0))
921 : END DO
922 : END IF
923 6 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
924 6 : 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 12 : CPASSERT(.NOT. do_kpoints)
929 12 : t0_state => wfi_get_snapshot(wf_history, wf_index=2)
930 12 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
931 12 : CPASSERT(ASSOCIATED(t0_state))
932 12 : CPASSERT(ASSOCIATED(t1_state))
933 12 : CPASSERT(ASSOCIATED(t0_state%wf))
934 12 : CPASSERT(ASSOCIATED(t1_state%wf))
935 12 : IF (wf_history%store_overlap) THEN
936 4 : CPASSERT(ASSOCIATED(t0_state%overlap))
937 4 : CPASSERT(ASSOCIATED(t1_state%overlap))
938 : END IF
939 12 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
940 12 : IF (nvec >= wf_history%memory_depth) THEN
941 12 : IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. (qs_env%scf_control%eps_scf_hist /= 0)) THEN
942 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
943 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
944 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
945 12 : ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
946 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
947 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
948 12 : ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
949 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
950 : END IF
951 : END IF
952 :
953 12 : 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 12 : CALL qs_rho_get(rho, rho_ao=rho_ao)
957 24 : DO ispin = 1, SIZE(mos)
958 12 : NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
959 12 : 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 12 : matrix_struct=matrix_struct)
962 : CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, &
963 12 : nrow_global=k, ncol_global=k)
964 12 : CALL cp_fm_create(csc, matrix_struct_new)
965 12 : CALL cp_fm_struct_release(matrix_struct_new)
966 :
967 12 : IF (use_overlap) THEN
968 4 : CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), mo_coeff, k)
969 4 : 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 8 : t1_state%wf(ispin), 0.0_dp, csc)
973 : END IF
974 12 : CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, mo_coeff)
975 12 : CALL cp_fm_release(csc)
976 12 : 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 12 : n_col=k)
980 : CALL calculate_density_matrix(mo_set=mos(ispin), &
981 48 : density_matrix=rho_ao(ispin)%matrix)
982 : END DO
983 12 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
984 12 : 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 376 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
989 376 : CPASSERT(nvec > 0)
990 376 : IF (nvec >= wf_history%memory_depth) THEN
991 178 : IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. (qs_env%scf_control%eps_scf_hist /= 0)) THEN
992 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
993 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
994 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
995 178 : ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
996 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
997 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
998 178 : ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
999 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1000 : END IF
1001 : END IF
1002 :
1003 376 : IF (do_kpoints) THEN
1004 4 : CALL wfi_extrapolate_ps_aspc_kp(wf_history, qs_env, nvec, io_unit, print_level)
1005 4 : my_orthogonal_wf = .TRUE.
1006 : ELSE
1007 372 : my_orthogonal_wf = .TRUE.
1008 822 : DO ispin = 1, SIZE(mos)
1009 450 : NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
1010 450 : 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 450 : matrix_struct=matrix_struct)
1013 450 : CALL cp_fm_create(fm_tmp, matrix_struct)
1014 : CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, &
1015 450 : nrow_global=k, ncol_global=k)
1016 450 : CALL cp_fm_create(csc, matrix_struct_new)
1017 450 : CALL cp_fm_struct_release(matrix_struct_new)
1018 : ! first the most recent
1019 450 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
1020 450 : CALL cp_fm_to_fm(t1_state%wf(ispin), mo_coeff)
1021 450 : alpha = nvec
1022 450 : CALL cp_fm_scale(alpha, mo_coeff)
1023 450 : CALL qs_rho_get(rho, rho_ao=rho_ao)
1024 962 : DO i = 2, nvec
1025 512 : t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1026 512 : IF (use_overlap) THEN
1027 474 : CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
1028 474 : 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 38 : t1_state%wf(ispin), 0.0_dp, csc)
1032 : END IF
1033 512 : CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
1034 512 : alpha = -1.0_dp*alpha*REAL(nvec - i + 1, dp)/REAL(i, dp)
1035 962 : CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp)
1036 : END DO
1037 :
1038 450 : CALL cp_fm_release(csc)
1039 450 : CALL cp_fm_release(fm_tmp)
1040 : CALL reorthogonalize_vectors(qs_env, &
1041 : v_matrix=mo_coeff, &
1042 450 : n_col=k)
1043 : CALL calculate_density_matrix(mo_set=mos(ispin), &
1044 1722 : density_matrix=rho_ao(ispin)%matrix)
1045 : END DO
1046 372 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1047 372 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
1048 : END IF
1049 :
1050 : CASE (wfi_aspc_nr)
1051 14888 : CALL cite_reference(Kolafa2004)
1052 14888 : CALL cite_reference(Kuhne2007)
1053 : ! figure out the actual number of vectors to use in the extrapolation:
1054 14888 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
1055 14888 : CPASSERT(nvec > 0)
1056 14888 : IF (nvec >= wf_history%memory_depth) THEN
1057 9658 : IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. &
1058 : (qs_env%scf_control%eps_scf_hist /= 0)) THEN
1059 18 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1060 18 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1061 18 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1062 9640 : ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
1063 62 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1064 62 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1065 9578 : ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
1066 8 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1067 : END IF
1068 : END IF
1069 :
1070 14888 : IF (do_kpoints) THEN
1071 420 : CALL wfi_extrapolate_ps_aspc_kp(wf_history, qs_env, nvec, io_unit, print_level)
1072 420 : my_orthogonal_wf = .TRUE.
1073 : ELSE
1074 14468 : my_orthogonal_wf = .TRUE.
1075 14468 : CALL qs_rho_get(rho, rho_ao=rho_ao)
1076 30155 : DO ispin = 1, SIZE(mos)
1077 15687 : NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
1078 15687 : 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 15687 : matrix_struct=matrix_struct)
1083 15687 : 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 15687 : ncol_global=k)
1088 15687 : CALL cp_fm_create(csc, matrix_struct_new, set_zero=.TRUE.)
1089 15687 : CALL cp_fm_struct_release(matrix_struct_new)
1090 : ! first the most recent
1091 : t1_state => wfi_get_snapshot(wf_history, &
1092 15687 : wf_index=1)
1093 15687 : CALL cp_fm_to_fm(t1_state%wf(ispin), mo_coeff)
1094 15687 : alpha = REAL(4*nvec - 2, KIND=dp)/REAL(nvec + 1, KIND=dp)
1095 15687 : 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 3138 : "Parameters for the always stable predictor-corrector (ASPC) method:", &
1098 3138 : "ASPC order: ", MAX(nvec - 2, 0), &
1099 6276 : "B(", 1, ") = ", alpha
1100 : END IF
1101 15687 : CALL cp_fm_scale(alpha, mo_coeff)
1102 :
1103 61279 : DO i = 2, nvec
1104 45592 : t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1105 45592 : IF (use_overlap) THEN
1106 34052 : CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
1107 34052 : 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 11540 : t1_state%wf(ispin), 0.0_dp, csc)
1111 : END IF
1112 45592 : 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 45592 : binomial(2*nvec, nvec - i)/binomial(2*nvec - 2, nvec - 1)
1115 45592 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1116 : WRITE (UNIT=io_unit, FMT="(T3,A2,I0,A4,F10.6)") &
1117 9414 : "B(", i, ") = ", alpha
1118 : END IF
1119 61279 : CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp)
1120 : END DO
1121 15687 : CALL cp_fm_release(csc)
1122 15687 : CALL cp_fm_release(fm_tmp)
1123 : CALL reorthogonalize_vectors(qs_env, &
1124 : v_matrix=mo_coeff, &
1125 15687 : n_col=k)
1126 : CALL calculate_density_matrix(mo_set=mos(ispin), &
1127 45842 : density_matrix=rho_ao(ispin)%matrix)
1128 : END DO
1129 14468 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1130 14468 : 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 142 : IF (do_kpoints) THEN
1135 6 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
1136 6 : CPASSERT(nvec > 0)
1137 6 : CALL wfi_extrapolate_gext_proj_kp(wf_history, qs_env, nvec, io_unit, print_level)
1138 6 : my_orthogonal_wf = .TRUE.
1139 : ELSE
1140 :
1141 : ! figure out the actual number of vectors to use in the extrapolation:
1142 136 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
1143 136 : IF (nvec >= wf_history%memory_depth) THEN
1144 88 : IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. &
1145 : (qs_env%scf_control%eps_scf_hist /= 0)) THEN
1146 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1147 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1148 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1149 88 : ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
1150 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1151 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1152 88 : ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
1153 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1154 : END IF
1155 : END IF
1156 136 : CPASSERT(nvec > 0)
1157 :
1158 : ! get the coefficients for the fitting
1159 408 : ALLOCATE (coeffs(nvec))
1160 136 : NULLIFY (matrix_s)
1161 136 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
1162 : CALL diff_fitting(wf_history, matrix_s(1)%matrix, coeffs, nvec, &
1163 136 : 1e-4_dp, io_unit, print_level)
1164 :
1165 136 : my_orthogonal_wf = .TRUE.
1166 136 : CALL qs_rho_get(rho, rho_ao=rho_ao)
1167 328 : DO ispin = 1, SIZE(mos)
1168 192 : NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
1169 192 : 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 192 : matrix_struct=matrix_struct)
1174 192 : 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 192 : ncol_global=k)
1179 192 : CALL cp_fm_create(csc, matrix_struct_new)
1180 192 : CALL cp_fm_struct_release(matrix_struct_new)
1181 :
1182 192 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
1183 :
1184 : ! do the linear combination of previous PSs
1185 192 : CALL cp_fm_set_all(mo_coeff, 0.0_dp)
1186 704 : DO i = 1, nvec
1187 512 : t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1188 512 : CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
1189 512 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
1190 512 : CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
1191 704 : CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, coeffs(i), fm_tmp)
1192 : END DO
1193 192 : CALL cp_fm_release(csc)
1194 192 : CALL cp_fm_release(fm_tmp)
1195 : CALL reorthogonalize_vectors(qs_env, &
1196 : v_matrix=mo_coeff, &
1197 192 : n_col=k)
1198 : CALL calculate_density_matrix(mo_set=mos(ispin), &
1199 712 : density_matrix=rho_ao(ispin)%matrix)
1200 : END DO
1201 136 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1202 136 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
1203 :
1204 136 : DEALLOCATE (coeffs)
1205 :
1206 : END IF
1207 :
1208 : CASE (wfi_gext_proj_qtr_nr)
1209 28 : IF (do_kpoints) THEN
1210 4 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
1211 4 : CPASSERT(nvec > 0)
1212 4 : CALL wfi_extrapolate_gext_proj_kp(wf_history, qs_env, nvec, io_unit, print_level)
1213 4 : my_orthogonal_wf = .TRUE.
1214 : ELSE
1215 :
1216 : ! figure out the actual number of vectors to use in the extrapolation:
1217 24 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
1218 24 : IF (nvec >= wf_history%memory_depth) THEN
1219 8 : IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. &
1220 : (qs_env%scf_control%eps_scf_hist /= 0)) THEN
1221 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1222 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1223 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1224 8 : ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
1225 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1226 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1227 8 : ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
1228 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1229 : END IF
1230 : END IF
1231 24 : CPASSERT(nvec > 0)
1232 :
1233 : ! get the coefficients for the fitting
1234 72 : ALLOCATE (coeffs(nvec))
1235 24 : NULLIFY (matrix_s)
1236 24 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
1237 : CALL tr_fitting(wf_history, matrix_s(1)%matrix, coeffs, nvec, &
1238 24 : 1e-4_dp, io_unit, print_level)
1239 :
1240 24 : my_orthogonal_wf = .TRUE.
1241 24 : CALL qs_rho_get(rho, rho_ao=rho_ao)
1242 48 : DO ispin = 1, SIZE(mos)
1243 24 : NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
1244 24 : 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 24 : matrix_struct=matrix_struct)
1249 24 : 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 24 : ncol_global=k)
1254 24 : CALL cp_fm_create(csc, matrix_struct_new)
1255 24 : CALL cp_fm_struct_release(matrix_struct_new)
1256 :
1257 24 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
1258 :
1259 : ! do the linear combination of previous PSs
1260 24 : CALL cp_fm_set_all(mo_coeff, 0.0_dp)
1261 104 : DO i = 1, nvec
1262 80 : t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1263 80 : CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
1264 80 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
1265 80 : CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
1266 104 : CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, coeffs(i), fm_tmp)
1267 : END DO
1268 24 : CALL cp_fm_release(csc)
1269 24 : CALL cp_fm_release(fm_tmp)
1270 : CALL reorthogonalize_vectors(qs_env, &
1271 : v_matrix=mo_coeff, &
1272 24 : n_col=k)
1273 : CALL calculate_density_matrix(mo_set=mos(ispin), &
1274 96 : density_matrix=rho_ao(ispin)%matrix)
1275 : END DO
1276 24 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1277 24 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
1278 :
1279 24 : DEALLOCATE (coeffs)
1280 :
1281 : END IF
1282 :
1283 : CASE default
1284 : CALL cp_abort(__LOCATION__, &
1285 : "Unknown interpolation method: "// &
1286 26671 : TRIM(ADJUSTL(cp_to_string(wf_history%interpolation_method_nr))))
1287 : END SELECT
1288 26671 : IF (PRESENT(orthogonal_wf)) orthogonal_wf = my_orthogonal_wf
1289 : CALL cp_print_key_finished_output(io_unit, logger, qs_env%input, &
1290 26671 : "DFT%SCF%PRINT%PROGRAM_RUN_INFO")
1291 26671 : CALL timestop(handle)
1292 26671 : 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 440 : 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 440 : 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 440 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: pbc_shift_cur, pbc_shift_src
1316 : INTEGER, DIMENSION(2) :: kp_range
1317 440 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
1318 440 : 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 440 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1323 440 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1324 440 : 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 440 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: csmat_cur
1328 440 : 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 440 : 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 440 : 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 440 : CALL timeset(routineN, handle)
1345 :
1346 440 : NULLIFY (kpoints, matrix_s_kp, scf_env, sab_nl, kp, &
1347 440 : mo_coeff, rmos, imos, wf_history, t1_state)
1348 :
1349 440 : 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 440 : sab_nl=sab_nl, cell_to_index=cell_to_index)
1353 440 : kplocal = kp_range(2) - kp_range(1) + 1
1354 :
1355 440 : IF (use_real_wfn) THEN
1356 0 : CALL timestop(handle)
1357 0 : RETURN
1358 : END IF
1359 :
1360 440 : wf_history => qs_env%wf_history
1361 440 : reload_snapshot_wf = .FALSE.
1362 440 : IF (PRESENT(load_snapshot_wf)) reload_snapshot_wf = load_snapshot_wf
1363 440 : IF (PRESENT(pbc_shift_ref)) THEN
1364 1302 : ALLOCATE (pbc_shift_src(3, SIZE(pbc_shift_ref, 2)))
1365 12034 : pbc_shift_src(:, :) = pbc_shift_ref(:, :)
1366 438 : use_pbc_phase_ref = .TRUE.
1367 : ELSE
1368 6 : use_pbc_phase_ref = .FALSE.
1369 6 : IF (ASSOCIATED(wf_history)) THEN
1370 6 : IF (wf_history%store_wf_kp .AND. wf_history%snapshot_count > 0) THEN
1371 4 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
1372 4 : CPASSERT(ASSOCIATED(t1_state%wf_kp))
1373 4 : CPASSERT(ASSOCIATED(t1_state%kp_pbc_shift))
1374 4 : reload_snapshot_wf = .TRUE.
1375 12 : ALLOCATE (pbc_shift_src(3, SIZE(t1_state%kp_pbc_shift, 2)))
1376 132 : pbc_shift_src(:, :) = t1_state%kp_pbc_shift(:, :)
1377 : use_pbc_phase_ref = .TRUE.
1378 : END IF
1379 : END IF
1380 : END IF
1381 438 : IF (use_pbc_phase_ref) CALL wfi_compute_kp_pbc_shift(qs_env, pbc_shift_cur)
1382 :
1383 440 : kp => kpoints%kp_env(1)%kpoint_env
1384 440 : nspin = SIZE(kp%mos, 2)
1385 440 : CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo, mo_coeff=mo_coeff)
1386 :
1387 440 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1388 : WRITE (UNIT=io_unit, FMT="(/,T2,A)") &
1389 0 : "Using previous wavefunctions as initial guess for k-points (with reorthogonalization)"
1390 : END IF
1391 :
1392 : ! Allocate dbcsr work matrices
1393 440 : ALLOCATE (rmatrix, cmatrix_db, tmpmat)
1394 440 : CALL dbcsr_create(rmatrix, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
1395 440 : CALL dbcsr_create(cmatrix_db, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
1396 440 : CALL dbcsr_create(tmpmat, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1397 440 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
1398 440 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix_db, sab_nl)
1399 :
1400 440 : CALL get_kpoint_info(kpoints, mpools=mpools_kp)
1401 440 : CALL mpools_get(mpools_kp, ao_ao_fm_pools=ao_ao_fm_pools_kp)
1402 440 : CALL fm_pool_create_fm(ao_ao_fm_pools_kp(1)%pool, fmlocal)
1403 440 : CALL cp_fm_get_info(fmlocal, matrix_struct=ao_ao_struct)
1404 :
1405 440 : CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct)
1406 440 : CALL cp_cfm_create(cfm_nao_nmo_work, mo_coeff%matrix_struct)
1407 :
1408 440 : NULLIFY (nmo_nmo_struct)
1409 : CALL cp_fm_struct_create(nmo_nmo_struct, template_fmstruct=mo_coeff%matrix_struct, &
1410 440 : nrow_global=nmo, ncol_global=nmo)
1411 440 : CALL cp_cfm_create(csc_cfm, nmo_nmo_struct)
1412 440 : CALL cp_fm_struct_release(nmo_nmo_struct)
1413 :
1414 440 : para_env => kpoints%blacs_env_all%para_env
1415 8764 : ALLOCATE (info(kplocal*nkp_groups, 2))
1416 :
1417 2458 : ALLOCATE (csmat_cur(kplocal))
1418 1578 : DO ikp = 1, kplocal
1419 1578 : 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 1578 : DO ikp = 1, kplocal
1425 3540 : DO igroup = 1, nkp_groups
1426 1962 : ik = kp_dist(1, igroup) + ikp - 1
1427 1962 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
1428 1962 : indx = indx + 1
1429 :
1430 1962 : CALL dbcsr_set(rmatrix, 0.0_dp)
1431 1962 : CALL dbcsr_set(cmatrix_db, 0.0_dp)
1432 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix_db, rsmat=matrix_s_kp, &
1433 1962 : ispin=1, xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1434 1962 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1435 1962 : CALL copy_dbcsr_to_fm(tmpmat, scf_env%scf_work1(1))
1436 1962 : CALL dbcsr_desymmetrize(cmatrix_db, tmpmat)
1437 1962 : CALL copy_dbcsr_to_fm(tmpmat, scf_env%scf_work1(2))
1438 :
1439 3100 : IF (my_kpgrp) THEN
1440 1138 : CALL cp_fm_start_copy_general(scf_env%scf_work1(1), fmlocal, para_env, info(indx, 1))
1441 1138 : CALL cp_fm_start_copy_general(scf_env%scf_work1(2), fmlocal, para_env, info(indx, 2))
1442 : ELSE
1443 824 : CALL cp_fm_start_copy_general(scf_env%scf_work1(1), fmdummy, para_env, info(indx, 1))
1444 824 : 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 1578 : DO ikp = 1, kplocal
1452 3540 : DO igroup = 1, nkp_groups
1453 1962 : ik = kp_dist(1, igroup) + ikp - 1
1454 1962 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
1455 824 : indx = indx + 1
1456 1138 : IF (my_kpgrp) THEN
1457 1138 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
1458 1138 : CALL cp_cfm_scale_and_add_fm(z_zero, csmat_cur(ikp), z_one, fmlocal)
1459 1138 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
1460 1138 : 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 2402 : DO indx = 1, kplocal*nkp_groups
1466 1962 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
1467 2402 : 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 1320 : ALLOCATE (eigenvalues(nmo))
1473 440 : eval_thresh = 1.0E-12_dp
1474 :
1475 1578 : DO ikp = 1, kplocal
1476 1138 : kp => kpoints%kp_env(ikp)%kpoint_env
1477 2472 : DO ispin = 1, nspin
1478 1334 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
1479 1334 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1480 1334 : IF (reload_snapshot_wf) THEN
1481 16 : CALL cp_fm_to_fm(t1_state%wf_kp(ikp, 1, ispin), rmos)
1482 16 : CALL cp_fm_to_fm(t1_state%wf_kp(ikp, 2, ispin), imos)
1483 : END IF
1484 1334 : IF (use_pbc_phase_ref) THEN
1485 1326 : ik = kp_range(1) + ikp - 1
1486 : CALL wfi_apply_kp_pbc_phase_fm(rmos, imos, pbc_shift_cur - pbc_shift_src, &
1487 25606 : xkp(1:3, ik), matrix_s_kp(1, 1)%matrix)
1488 : END IF
1489 1334 : CALL cp_fm_to_cfm(rmos, imos, cmos_new)
1490 :
1491 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nao, z_one, &
1492 1334 : csmat_cur(ikp), cmos_new, z_zero, cfm_nao_nmo_work)
1493 : CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, &
1494 1334 : cmos_new, cfm_nao_nmo_work, z_zero, csc_cfm)
1495 :
1496 1334 : CALL cp_cfm_cholesky_decompose(csc_cfm, info_out=chol_info)
1497 1334 : IF (chol_info == 0) THEN
1498 1334 : CALL cp_cfm_triangular_multiply(csc_cfm, cmos_new, side='R', invert_tr=.TRUE., uplo_tr='U')
1499 : ELSE
1500 0 : CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, cmos_new, cfm_nao_nmo_work, z_zero, csc_cfm)
1501 0 : CALL cp_cfm_create(cfm_evecs, csc_cfm%matrix_struct)
1502 0 : CALL cp_cfm_create(cfm_mhalf, csc_cfm%matrix_struct)
1503 0 : CALL cp_cfm_heevd(csc_cfm, cfm_evecs, eigenvalues)
1504 0 : CALL cp_cfm_to_cfm(cfm_evecs, cfm_mhalf)
1505 0 : ALLOCATE (col_scaling(nmo))
1506 0 : DO j = 1, nmo
1507 0 : IF (eigenvalues(j) > eval_thresh) THEN
1508 0 : col_scaling(j) = CMPLX(1.0_dp/SQRT(eigenvalues(j)), 0.0_dp, KIND=dp)
1509 : ELSE
1510 0 : col_scaling(j) = z_zero
1511 : END IF
1512 : END DO
1513 0 : CALL cp_cfm_column_scale(cfm_mhalf, col_scaling)
1514 0 : DEALLOCATE (col_scaling)
1515 0 : CALL cp_cfm_gemm('N', 'C', nmo, nmo, nmo, z_one, cfm_mhalf, cfm_evecs, z_zero, csc_cfm)
1516 0 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nmo, z_one, cmos_new, csc_cfm, z_zero, cfm_nao_nmo_work)
1517 0 : CALL cp_cfm_to_cfm(cfm_nao_nmo_work, cmos_new)
1518 0 : CALL cp_cfm_release(cfm_evecs)
1519 0 : CALL cp_cfm_release(cfm_mhalf)
1520 : END IF
1521 3806 : 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 1578 : kp%mos_prefilled = .TRUE.
1529 : END DO
1530 440 : DEALLOCATE (eigenvalues)
1531 :
1532 : ! Phase C: Rebuild Density Matrix P(R)
1533 440 : CALL qs_kpoint_state_commit(qs_env, update_occupations=.TRUE.)
1534 :
1535 : ! Cleanup
1536 1578 : DO ikp = 1, kplocal
1537 1578 : CALL cp_cfm_release(csmat_cur(ikp))
1538 : END DO
1539 440 : DEALLOCATE (csmat_cur)
1540 4364 : DEALLOCATE (info)
1541 440 : CALL cp_cfm_release(cmos_new)
1542 440 : CALL cp_cfm_release(cfm_nao_nmo_work)
1543 440 : CALL cp_cfm_release(csc_cfm)
1544 440 : CALL fm_pool_give_back_fm(ao_ao_fm_pools_kp(1)%pool, fmlocal)
1545 440 : CALL dbcsr_deallocate_matrix(rmatrix)
1546 440 : CALL dbcsr_deallocate_matrix(cmatrix_db)
1547 440 : CALL dbcsr_deallocate_matrix(tmpmat)
1548 440 : IF (ALLOCATED(pbc_shift_cur)) DEALLOCATE (pbc_shift_cur)
1549 440 : IF (ALLOCATED(pbc_shift_src)) DEALLOCATE (pbc_shift_src)
1550 :
1551 440 : CALL timestop(handle)
1552 2640 : 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 2926 : 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 2926 : CPASSERT(ASSOCIATED(snapshot))
1570 2926 : CPASSERT(ASSOCIATED(cell))
1571 2926 : CPASSERT(ASSOCIATED(particle_set))
1572 :
1573 2926 : natom = SIZE(particle_set)
1574 2926 : IF (ASSOCIATED(snapshot%kp_pbc_shift)) THEN
1575 300 : DEALLOCATE (snapshot%kp_pbc_shift)
1576 : END IF
1577 8778 : ALLOCATE (snapshot%kp_pbc_shift(3, natom))
1578 15804 : DO iatom = 1, natom
1579 12878 : r_pbc(1:3) = pbc(particle_set(iatom)%r(1:3), cell)
1580 12878 : CALL real_to_scaled(frac_raw, particle_set(iatom)%r(1:3), cell)
1581 12878 : CALL real_to_scaled(frac_pbc, r_pbc(1:3), cell)
1582 54438 : snapshot%kp_pbc_shift(1:3, iatom) = NINT(frac_pbc(1:3) - frac_raw(1:3))
1583 : END DO
1584 2926 : 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 438 : 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 438 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1599 :
1600 438 : NULLIFY (cell, particle_set)
1601 438 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1602 438 : CPASSERT(ASSOCIATED(cell))
1603 438 : CPASSERT(ASSOCIATED(particle_set))
1604 :
1605 438 : natom = SIZE(particle_set)
1606 1314 : ALLOCATE (pbc_shift(3, natom))
1607 3370 : DO iatom = 1, natom
1608 2932 : r_pbc(1:3) = pbc(particle_set(iatom)%r(1:3), cell)
1609 2932 : CALL real_to_scaled(frac_raw, particle_set(iatom)%r(1:3), cell)
1610 2932 : CALL real_to_scaled(frac_pbc, r_pbc(1:3), cell)
1611 12166 : pbc_shift(1:3, iatom) = NINT(frac_pbc(1:3) - frac_raw(1:3))
1612 : END DO
1613 438 : 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 1326 : 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 1326 : INTEGER, DIMENSION(:), POINTER :: row_blk_size
1633 : REAL(KIND=dp) :: ci, cr, i_old, r_old, theta
1634 1326 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: iblock, rblock
1635 :
1636 0 : CPASSERT(ASSOCIATED(rmos))
1637 1326 : CPASSERT(ASSOCIATED(imos))
1638 1326 : CPASSERT(ASSOCIATED(matrix_template))
1639 :
1640 1326 : natom = SIZE(pbc_shift_delta, 2)
1641 1326 : CALL cp_fm_get_info(rmos, ncol_global=nmo)
1642 1326 : NULLIFY (row_blk_size)
1643 1326 : CALL dbcsr_get_info(matrix_template, row_blk_size=row_blk_size)
1644 1326 : CPASSERT(SIZE(row_blk_size) >= natom)
1645 :
1646 1326 : row_start = 1
1647 7396 : DO iatom = 1, natom
1648 6070 : nrow = row_blk_size(iatom)
1649 24056 : IF (ANY(pbc_shift_delta(1:3, iatom) /= 0)) THEN
1650 528 : theta = twopi*SUM(xk(1:3)*REAL(pbc_shift_delta(1:3, iatom), KIND=dp))
1651 132 : cr = COS(theta)
1652 132 : ci = SIN(theta)
1653 792 : ALLOCATE (rblock(nrow, nmo), iblock(nrow, nmo))
1654 132 : CALL cp_fm_get_submatrix(rmos, rblock, row_start, 1, nrow, nmo)
1655 132 : CALL cp_fm_get_submatrix(imos, iblock, row_start, 1, nrow, nmo)
1656 4956 : DO icol = 1, nmo
1657 36612 : DO irow = 1, nrow
1658 31656 : r_old = rblock(irow, icol)
1659 31656 : i_old = iblock(irow, icol)
1660 31656 : rblock(irow, icol) = cr*r_old - ci*i_old
1661 36480 : iblock(irow, icol) = ci*r_old + cr*i_old
1662 : END DO
1663 : END DO
1664 132 : CALL cp_fm_set_submatrix(rmos, rblock, row_start, 1, nrow, nmo)
1665 132 : CALL cp_fm_set_submatrix(imos, iblock, row_start, 1, nrow, nmo)
1666 132 : DEALLOCATE (rblock, iblock)
1667 : END IF
1668 7396 : row_start = row_start + nrow
1669 : END DO
1670 2652 : 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 5348 : 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 5348 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: zblock
1688 : INTEGER :: iatom, natom, nmo, nrow, row_start
1689 5348 : INTEGER, DIMENSION(:), POINTER :: row_blk_size
1690 : REAL(KIND=dp) :: theta
1691 :
1692 0 : CPASSERT(ASSOCIATED(matrix_template))
1693 :
1694 5348 : natom = SIZE(pbc_shift_delta, 2)
1695 5348 : CALL cp_cfm_get_info(cmos, ncol_global=nmo)
1696 5348 : NULLIFY (row_blk_size)
1697 5348 : CALL dbcsr_get_info(matrix_template, row_blk_size=row_blk_size)
1698 5348 : CPASSERT(SIZE(row_blk_size) >= natom)
1699 :
1700 5348 : row_start = 1
1701 37252 : DO iatom = 1, natom
1702 31904 : nrow = row_blk_size(iatom)
1703 126784 : IF (ANY(pbc_shift_delta(1:3, iatom) /= 0)) THEN
1704 1856 : theta = twopi*SUM(xk(1:3)*REAL(pbc_shift_delta(1:3, iatom), KIND=dp))
1705 464 : phase = CMPLX(COS(theta), SIN(theta), KIND=dp)
1706 1856 : ALLOCATE (zblock(nrow, nmo))
1707 464 : CALL cp_cfm_get_submatrix(cmos, zblock, row_start, 1, nrow, nmo)
1708 197224 : zblock = phase*zblock
1709 464 : CALL cp_cfm_set_submatrix(cmos, zblock, row_start, 1, nrow, nmo)
1710 464 : DEALLOCATE (zblock)
1711 : END IF
1712 37252 : row_start = row_start + nrow
1713 : END DO
1714 10696 : 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 848 : 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 424 : 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 424 : 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 424 : method_nr = wf_history%interpolation_method_nr
1750 :
1751 424 : CALL timeset(routineN, handle)
1752 424 : NULLIFY (kpoints, kp, mo_coeff, rmos, imos, t0_state, t1_state, nmo_nmo_struct, &
1753 424 : matrix_s_kp, xkp)
1754 :
1755 424 : CALL get_qs_env(qs_env, kpoints=kpoints, matrix_s_kp=matrix_s_kp)
1756 424 : CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn, kp_range=kp_range, xkp=xkp)
1757 424 : kplocal = kp_range(2) - kp_range(1) + 1
1758 :
1759 424 : IF (use_real_wfn) THEN
1760 0 : IF (method_nr == wfi_aspc_nr) THEN
1761 : CALL cp_warn(__LOCATION__, "ASPC with k-points requires complex wavefunctions; "// &
1762 0 : "falling back to USE_PREV_WF.")
1763 : ELSE
1764 : CALL cp_warn(__LOCATION__, "PS with k-points requires complex wavefunctions; "// &
1765 0 : "falling back to USE_PREV_WF.")
1766 : END IF
1767 0 : CALL wfi_use_prev_wf_kp(qs_env, io_unit, print_level)
1768 0 : CALL timestop(handle)
1769 0 : RETURN
1770 : END IF
1771 :
1772 424 : kp => kpoints%kp_env(1)%kpoint_env
1773 424 : nspin = SIZE(kp%mos, 2)
1774 424 : CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo, mo_coeff=mo_coeff)
1775 :
1776 424 : IF (method_nr == wfi_aspc_nr) THEN
1777 420 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1778 : WRITE (UNIT=io_unit, FMT="(/,T2,A,/,T3,A,I0)") &
1779 40 : "Parameters for the always stable predictor-corrector (ASPC) method:", &
1780 80 : "ASPC order: ", MAX(nvec - 2, 0)
1781 : END IF
1782 : END IF
1783 :
1784 16 : IF (method_nr == wfi_aspc_nr) THEN
1785 420 : CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct, set_zero=.TRUE.)
1786 420 : CALL cp_cfm_create(cmos_1, mo_coeff%matrix_struct, set_zero=.TRUE.)
1787 420 : CALL cp_cfm_create(cmos_i, mo_coeff%matrix_struct, set_zero=.TRUE.)
1788 420 : CALL cp_cfm_create(cfm_nao_nmo_work, mo_coeff%matrix_struct, set_zero=.TRUE.)
1789 : ELSE
1790 4 : CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct)
1791 4 : CALL cp_cfm_create(cmos_1, mo_coeff%matrix_struct)
1792 4 : CALL cp_cfm_create(cmos_i, mo_coeff%matrix_struct)
1793 4 : 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 424 : nrow_global=nmo, ncol_global=nmo)
1798 424 : IF (method_nr == wfi_aspc_nr) THEN
1799 420 : CALL cp_cfm_create(csc_cfm, nmo_nmo_struct, set_zero=.TRUE.)
1800 : ELSE
1801 4 : CALL cp_cfm_create(csc_cfm, nmo_nmo_struct)
1802 : END IF
1803 424 : CALL cp_fm_struct_release(nmo_nmo_struct)
1804 :
1805 : ! Phase 1: Initialize C_new(k) = B(1) * C_1(k)
1806 424 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
1807 424 : IF (method_nr == wfi_aspc_nr) THEN
1808 420 : alpha_coeff = REAL(4*nvec - 2, KIND=dp)/REAL(nvec + 1, KIND=dp)
1809 420 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1810 40 : WRITE (UNIT=io_unit, FMT="(T3,A2,I0,A4,F10.6)") "B(", 1, ") = ", alpha_coeff
1811 : END IF
1812 : ELSE
1813 4 : alpha_coeff = nvec
1814 : END IF
1815 :
1816 1504 : DO ikp = 1, kplocal
1817 1080 : kp => kpoints%kp_env(ikp)%kpoint_env
1818 2780 : DO ispin = 1, nspin
1819 1276 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
1820 1276 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1821 1276 : CALL cp_fm_to_fm(t1_state%wf_kp(ikp, 1, ispin), rmos)
1822 1276 : CALL cp_fm_to_fm(t1_state%wf_kp(ikp, 2, ispin), imos)
1823 1276 : CALL cp_fm_scale(alpha_coeff, rmos)
1824 2356 : 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 1736 : DO i = 2, nvec
1830 1312 : t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1831 1312 : IF (method_nr == wfi_aspc_nr) THEN
1832 : alpha_coeff = (-1.0_dp)**(i + 1)*REAL(i, KIND=dp)* &
1833 1310 : binomial(2*nvec, nvec - i)/binomial(2*nvec - 2, nvec - 1)
1834 1310 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1835 71 : WRITE (UNIT=io_unit, FMT="(T3,A2,I0,A4,F10.6)") "B(", i, ") = ", alpha_coeff
1836 : END IF
1837 : ELSE
1838 2 : alpha_coeff = -1.0_dp*alpha_coeff*REAL(nvec - i + 1, dp)/REAL(i, dp)
1839 : END IF
1840 :
1841 4264 : DO ikp = 1, kplocal
1842 2528 : kp => kpoints%kp_env(ikp)%kpoint_env
1843 6464 : DO ispin = 1, nspin
1844 2624 : ik = kp_range(1) + ikp - 1
1845 2624 : CALL cp_fm_to_cfm(t1_state%wf_kp(ikp, 1, ispin), t1_state%wf_kp(ikp, 2, ispin), cmos_1)
1846 2624 : 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 64888 : 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 2624 : 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 2624 : cmos_i, cfm_nao_nmo_work, z_zero, csc_cfm)
1858 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nmo, z_one, &
1859 2624 : 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 64888 : xkp(1:3, ik), matrix_s_kp(1, 1)%matrix)
1866 :
1867 2624 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
1868 2624 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1869 2624 : CALL cp_fm_to_cfm(rmos, imos, cmos_new)
1870 2624 : CALL cp_cfm_scale_and_add(z_one, cmos_new, CMPLX(alpha_coeff, 0.0_dp, KIND=dp), cfm_nao_nmo_work)
1871 5152 : CALL cp_cfm_to_fm(cmos_new, rmos, imos)
1872 : END DO
1873 : END DO
1874 : END DO
1875 :
1876 424 : CALL cp_cfm_release(cmos_new)
1877 424 : CALL cp_cfm_release(cmos_1)
1878 424 : CALL cp_cfm_release(cmos_i)
1879 424 : CALL cp_cfm_release(cfm_nao_nmo_work)
1880 424 : 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 424 : load_snapshot_wf=.FALSE.)
1888 :
1889 424 : CALL timestop(handle)
1890 :
1891 424 : 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 10 : 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 10 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
1915 10 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1916 : LOGICAL :: my_kpgrp, use_real_wfn
1917 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coeffs, weight_kp
1918 10 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
1919 10 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1920 10 : 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 10 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: csmat_cur
1924 10 : 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 10 : 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 10 : 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 10 : method_nr = wf_history%interpolation_method_nr
1940 :
1941 10 : CALL timeset(routineN, handle)
1942 10 : NULLIFY (ao_ao_struct, cell_to_index, cmatrix_db, imos, kp, kpoints, matrix_s_kp, &
1943 10 : mo_coeff, mpools_kp, para_env, para_env_inter_kp, rmatrix, rmos, sab_nl, &
1944 10 : scf_env, t0_state, t1_state, tmpmat, xkp, kp_dist, wkp, nmo_nmo_struct, &
1945 10 : ao_ao_fm_pools_kp)
1946 :
1947 10 : 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 10 : mpools=mpools_kp, para_env_inter_kp=para_env_inter_kp)
1952 10 : kplocal = kp_range(2) - kp_range(1) + 1
1953 :
1954 10 : IF (use_real_wfn) THEN
1955 : CALL cp_warn(__LOCATION__, "GExt with k-points requires complex wavefunctions; "// &
1956 0 : "falling back to USE_PREV_WF.")
1957 0 : CALL wfi_use_prev_wf_kp(qs_env, io_unit, print_level)
1958 0 : CALL timestop(handle)
1959 0 : RETURN
1960 : END IF
1961 :
1962 10 : IF (nvec >= wf_history%memory_depth) THEN
1963 0 : IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. &
1964 : (qs_env%scf_control%eps_scf_hist /= 0)) THEN
1965 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1966 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1967 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1968 0 : ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
1969 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1970 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1971 0 : ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
1972 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1973 : END IF
1974 : END IF
1975 :
1976 10 : kp => kpoints%kp_env(1)%kpoint_env
1977 10 : nspin = SIZE(kp%mos, 2)
1978 10 : 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 10 : ALLOCATE (rmatrix, cmatrix_db, tmpmat)
1983 10 : CALL dbcsr_create(rmatrix, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
1984 10 : CALL dbcsr_create(cmatrix_db, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
1985 10 : CALL dbcsr_create(tmpmat, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1986 10 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
1987 10 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix_db, sab_nl)
1988 :
1989 10 : CALL mpools_get(mpools_kp, ao_ao_fm_pools=ao_ao_fm_pools_kp)
1990 10 : CALL fm_pool_create_fm(ao_ao_fm_pools_kp(1)%pool, fmlocal)
1991 10 : CALL cp_fm_get_info(fmlocal, matrix_struct=ao_ao_struct)
1992 :
1993 10 : para_env => kpoints%blacs_env_all%para_env
1994 178 : ALLOCATE (info(kplocal*nkp_groups, 2))
1995 84 : ALLOCATE (csmat_cur(kplocal), weight_kp(kplocal))
1996 44 : DO ikp = 1, kplocal
1997 34 : CALL cp_cfm_create(csmat_cur(ikp), ao_ao_struct)
1998 44 : weight_kp(ikp) = wkp(kp_range(1) + ikp - 1)
1999 : END DO
2000 :
2001 : indx = 0
2002 44 : DO ikp = 1, kplocal
2003 78 : DO igroup = 1, nkp_groups
2004 34 : ik = kp_dist(1, igroup) + ikp - 1
2005 34 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
2006 34 : indx = indx + 1
2007 :
2008 34 : CALL dbcsr_set(rmatrix, 0.0_dp)
2009 34 : CALL dbcsr_set(cmatrix_db, 0.0_dp)
2010 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix_db, rsmat=matrix_s_kp, &
2011 34 : ispin=1, xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
2012 34 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
2013 34 : CALL copy_dbcsr_to_fm(tmpmat, scf_env%scf_work1(1))
2014 34 : CALL dbcsr_desymmetrize(cmatrix_db, tmpmat)
2015 34 : CALL copy_dbcsr_to_fm(tmpmat, scf_env%scf_work1(2))
2016 :
2017 68 : IF (my_kpgrp) THEN
2018 34 : CALL cp_fm_start_copy_general(scf_env%scf_work1(1), fmlocal, para_env, info(indx, 1))
2019 34 : CALL cp_fm_start_copy_general(scf_env%scf_work1(2), fmlocal, para_env, info(indx, 2))
2020 : ELSE
2021 0 : CALL cp_fm_start_copy_general(scf_env%scf_work1(1), fmdummy, para_env, info(indx, 1))
2022 0 : 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 44 : DO ikp = 1, kplocal
2029 78 : DO igroup = 1, nkp_groups
2030 34 : ik = kp_dist(1, igroup) + ikp - 1
2031 34 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
2032 0 : indx = indx + 1
2033 34 : IF (my_kpgrp) THEN
2034 34 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
2035 34 : CALL cp_cfm_scale_and_add_fm(z_zero, csmat_cur(ikp), z_one, fmlocal)
2036 34 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
2037 34 : 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 44 : DO indx = 1, kplocal*nkp_groups
2043 34 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
2044 44 : CALL cp_fm_cleanup_copy_general(info(indx, 2))
2045 : END DO
2046 :
2047 30 : ALLOCATE (coeffs(nvec))
2048 10 : 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 6 : 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 4 : 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 10 : CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct)
2060 10 : CALL cp_cfm_create(cmos_1, mo_coeff%matrix_struct)
2061 10 : CALL cp_cfm_create(cmos_i, mo_coeff%matrix_struct)
2062 10 : 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 10 : nrow_global=nmo, ncol_global=nmo)
2065 10 : CALL cp_cfm_create(csc_cfm, nmo_nmo_struct)
2066 10 : CALL cp_fm_struct_release(nmo_nmo_struct)
2067 :
2068 10 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
2069 44 : DO ikp = 1, kplocal
2070 34 : kp => kpoints%kp_env(ikp)%kpoint_env
2071 78 : DO ispin = 1, nspin
2072 34 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
2073 34 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
2074 34 : CALL cp_fm_set_all(rmos, 0.0_dp)
2075 68 : CALL cp_fm_set_all(imos, 0.0_dp)
2076 : END DO
2077 : END DO
2078 :
2079 24 : DO i = 1, nvec
2080 14 : t0_state => wfi_get_snapshot(wf_history, wf_index=i)
2081 74 : DO ikp = 1, kplocal
2082 50 : kp => kpoints%kp_env(ikp)%kpoint_env
2083 50 : ik = kp_range(1) + ikp - 1
2084 114 : DO ispin = 1, nspin
2085 50 : CALL cp_fm_to_cfm(t1_state%wf_kp(ikp, 1, ispin), t1_state%wf_kp(ikp, 2, ispin), cmos_1)
2086 50 : 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 1594 : xkp(1:3, ik), matrix_s_kp(1, 1)%matrix)
2090 :
2091 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nao, z_one, &
2092 50 : 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 50 : cmos_i, cfm_nao_nmo_work, z_zero, csc_cfm)
2095 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nmo, z_one, &
2096 50 : 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 1594 : xkp(1:3, ik), matrix_s_kp(1, 1)%matrix)
2100 :
2101 50 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
2102 50 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
2103 50 : CALL cp_fm_to_cfm(rmos, imos, cmos_new)
2104 50 : CALL cp_cfm_scale_and_add(z_one, cmos_new, CMPLX(coeffs(i), 0.0_dp, KIND=dp), cfm_nao_nmo_work)
2105 100 : CALL cp_cfm_to_fm(cmos_new, rmos, imos)
2106 : END DO
2107 : END DO
2108 : END DO
2109 :
2110 10 : CALL cp_cfm_release(cmos_new)
2111 10 : CALL cp_cfm_release(cmos_1)
2112 10 : CALL cp_cfm_release(cmos_i)
2113 10 : CALL cp_cfm_release(cfm_nao_nmo_work)
2114 10 : 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 10 : load_snapshot_wf=.FALSE.)
2118 :
2119 44 : DO ikp = 1, kplocal
2120 44 : CALL cp_cfm_release(csmat_cur(ikp))
2121 : END DO
2122 78 : DEALLOCATE (csmat_cur, coeffs, weight_kp, info)
2123 10 : CALL fm_pool_give_back_fm(ao_ao_fm_pools_kp(1)%pool, fmlocal)
2124 10 : CALL dbcsr_deallocate_matrix(rmatrix)
2125 10 : CALL dbcsr_deallocate_matrix(cmatrix_db)
2126 10 : CALL dbcsr_deallocate_matrix(tmpmat)
2127 :
2128 10 : CALL timestop(handle)
2129 :
2130 60 : 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 11225 : 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 11225 : IF (nvec >= qs_env%wf_history%memory_depth) THEN
2146 1693 : IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. (qs_env%scf_control%eps_scf_hist /= 0)) THEN
2147 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
2148 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
2149 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
2150 1693 : ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
2151 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
2152 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
2153 1693 : ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
2154 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
2155 0 : qs_env%scf_control%outer_scf%eps_scf = qs_env%scf_control%eps_scf_hist
2156 : END IF
2157 : END IF
2158 :
2159 11225 : 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 26675 : 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 26675 : CPASSERT(ASSOCIATED(wf_history))
2176 26675 : CPASSERT(wf_history%ref_count > 0)
2177 26675 : CPASSERT(ASSOCIATED(qs_env))
2178 :
2179 26675 : wf_history%snapshot_count = wf_history%snapshot_count + 1
2180 26675 : IF (wf_history%memory_depth > 0) THEN
2181 : wf_history%last_state_index = MODULO(wf_history%snapshot_count, &
2182 25480 : 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 25480 : qs_env=qs_env, dt=dt)
2186 : END IF
2187 26675 : 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 33778 : 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 16889 : 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 16889 : NULLIFY (scf_env, scf_control, maxao_maxmo_fm_pool, matrix_s, mpools, dft_control)
2217 16889 : CALL timeset(routineN, handle)
2218 :
2219 16889 : CPASSERT(ASSOCIATED(qs_env))
2220 :
2221 16889 : CALL cp_fm_get_info(v_matrix, ncol_global=my_n_col)
2222 16889 : 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 16889 : dft_control=dft_control)
2228 16889 : CALL mpools_get(mpools, maxao_maxmo_fm_pool=maxao_maxmo_fm_pool)
2229 16889 : IF (ASSOCIATED(scf_env)) THEN
2230 : ortho_contains_cholesky = (scf_env%method /= ot_method_nr) .AND. &
2231 : (scf_env%cholesky_method > 0) .AND. &
2232 16889 : ASSOCIATED(scf_env%ortho)
2233 : ELSE
2234 : ortho_contains_cholesky = .FALSE.
2235 : END IF
2236 :
2237 16889 : CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
2238 16889 : smearing_is_used = .FALSE.
2239 16889 : IF (dft_control%smear) THEN
2240 1906 : smearing_is_used = .TRUE.
2241 : END IF
2242 :
2243 16889 : IF (has_unit_metric) THEN
2244 3410 : CALL make_basis_simple(v_matrix, my_n_col)
2245 13479 : ELSE IF (smearing_is_used) THEN
2246 : CALL make_basis_lowdin(vmatrix=v_matrix, ncol=my_n_col, &
2247 1906 : matrix_s=matrix_s(1)%matrix)
2248 11573 : ELSE IF (ortho_contains_cholesky) THEN
2249 : CALL make_basis_cholesky(vmatrix=v_matrix, ncol=my_n_col, &
2250 8176 : ortho=scf_env%ortho)
2251 : ELSE
2252 3397 : CALL make_basis_sm(v_matrix, my_n_col, matrix_s(1)%matrix)
2253 : END IF
2254 16889 : CALL timestop(handle)
2255 16889 : 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 0 : 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 0 : NULLIFY (dft_control, wf_history)
2276 :
2277 0 : CALL timeset(routineN, handle)
2278 0 : logger => cp_get_default_logger()
2279 0 : 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 0 : extension=".scfLog")
2282 :
2283 0 : CPASSERT(ASSOCIATED(qs_env))
2284 0 : CPASSERT(ASSOCIATED(qs_env%wf_history))
2285 0 : CPASSERT(qs_env%wf_history%ref_count > 0)
2286 0 : CALL get_qs_env(qs_env, dft_control=dft_control)
2287 :
2288 0 : 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 0 : IF (qs_env%wf_history%snapshot_count >= 2) THEN
2296 0 : IF (debug_this_module .AND. io_unit > 0) THEN
2297 0 : 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 0 : has_unit_metric=qs_env%has_unit_metric)
2303 : CALL set_qs_env(qs_env=qs_env, &
2304 0 : wf_history=wf_history)
2305 0 : CALL wfi_release(wf_history)
2306 0 : CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
2307 : END IF
2308 : CASE DEFAULT
2309 0 : CPABORT("Unknown extrapolation method.")
2310 : END SELECT
2311 0 : CALL timestop(handle)
2312 :
2313 0 : 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 142 : SUBROUTINE diff_fitting(wf_history, current_overlap, coeffs, nvec, eps, io_unit, print_level, &
2335 142 : 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 142 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: b
2351 142 : 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 142 : IF (nvec <= 0) THEN
2358 0 : CPABORT("Not enough vectors to do the fitting")
2359 142 : ELSE IF (nvec == 1) THEN
2360 24 : coeffs(1) = 1.0_dp
2361 76 : RETURN
2362 : END IF
2363 :
2364 118 : IF (PRESENT(current_overlap_kp)) THEN
2365 12 : ALLOCATE (A(nvec - 1, nvec - 1), b(nvec - 1))
2366 2 : A = 0.0_dp
2367 2 : b = 0.0_dp
2368 :
2369 2 : ref_state => wfi_get_snapshot(wf_history, wf_index=1)
2370 2 : CALL cp_cfm_create(target_diff_cfm, current_overlap_kp(1)%matrix_struct)
2371 2 : CALL cp_cfm_create(tmp_i_cfm, current_overlap_kp(1)%matrix_struct)
2372 2 : CALL cp_cfm_create(tmp_j_cfm, current_overlap_kp(1)%matrix_struct)
2373 2 : CALL cp_cfm_create(tmp_conj_cfm, current_overlap_kp(1)%matrix_struct)
2374 :
2375 10 : DO ikp = 1, SIZE(current_overlap_kp)
2376 8 : weight = 1.0_dp
2377 8 : IF (PRESENT(kpoint_weights)) weight = kpoint_weights(ikp)
2378 :
2379 8 : 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 8 : ref_state%overlap_cfm_kp(ikp))
2382 18 : DO i = 2, nvec
2383 8 : state => wfi_get_snapshot(wf_history, wf_index=i)
2384 8 : 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 8 : ref_state%overlap_cfm_kp(ikp))
2387 8 : CALL cp_cfm_to_cfm(tmp_i_cfm, tmp_conj_cfm)
2388 264 : DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
2389 4360 : DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
2390 : tmp_conj_cfm%local_data(irow_local, icol_local) = &
2391 4352 : CONJG(tmp_conj_cfm%local_data(irow_local, icol_local))
2392 : END DO
2393 : END DO
2394 8 : CALL cp_cfm_trace(tmp_conj_cfm, target_diff_cfm, ztrace)
2395 8 : b(i - 1) = b(i - 1) + weight*REAL(ztrace, KIND=dp)
2396 :
2397 24 : DO j = 2, i
2398 8 : state => wfi_get_snapshot(wf_history, wf_index=j)
2399 8 : 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 8 : ref_state%overlap_cfm_kp(ikp))
2402 8 : CALL cp_cfm_trace(tmp_conj_cfm, tmp_j_cfm, ztrace)
2403 16 : 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 4 : DO i = 2, nvec
2409 6 : DO j = 2, i
2410 4 : A(i - 1, j - 1) = A(j - 1, i - 1)
2411 : END DO
2412 : END DO
2413 :
2414 2 : IF (PRESENT(para_env_inter_kp)) THEN
2415 2 : IF (ASSOCIATED(para_env_inter_kp)) THEN
2416 2 : CALL para_env_inter_kp%sum(A)
2417 2 : CALL para_env_inter_kp%sum(b)
2418 : END IF
2419 : END IF
2420 :
2421 4 : DO i = 1, nvec - 1
2422 4 : A(i, i) = A(i, i) + eps**2
2423 : END DO
2424 :
2425 2 : CALL dposv('u', nvec - 1, 1, A, nvec - 1, b, nvec - 1, info)
2426 2 : IF (info /= 0) THEN
2427 0 : CPABORT("DPOSV failed.")
2428 : END IF
2429 :
2430 4 : coeffs(1) = 1.0_dp - SUM(b)
2431 4 : coeffs(2:nvec) = b(:)
2432 :
2433 2 : IF (print_level > low_print_level) THEN
2434 2 : error = 0.0_dp
2435 2 : norm_ref = 0.0_dp
2436 10 : DO ikp = 1, SIZE(current_overlap_kp)
2437 8 : weight = 1.0_dp
2438 8 : IF (PRESENT(kpoint_weights)) weight = kpoint_weights(ikp)
2439 8 : CALL cp_cfm_to_cfm(current_overlap_kp(ikp), tmp_i_cfm)
2440 24 : DO i = 1, nvec
2441 16 : 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 24 : state%overlap_cfm_kp(ikp))
2444 : END DO
2445 8 : CALL cp_cfm_to_cfm(tmp_i_cfm, tmp_conj_cfm)
2446 264 : DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
2447 4360 : DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
2448 : tmp_conj_cfm%local_data(irow_local, icol_local) = &
2449 4352 : CONJG(tmp_conj_cfm%local_data(irow_local, icol_local))
2450 : END DO
2451 : END DO
2452 8 : CALL cp_cfm_trace(tmp_conj_cfm, tmp_i_cfm, ztrace)
2453 8 : error = error + weight*REAL(ztrace, KIND=dp)
2454 8 : CALL cp_cfm_to_cfm(ref_state%overlap_cfm_kp(ikp), tmp_conj_cfm)
2455 264 : DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
2456 4360 : DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
2457 : tmp_conj_cfm%local_data(irow_local, icol_local) = &
2458 4352 : CONJG(tmp_conj_cfm%local_data(irow_local, icol_local))
2459 : END DO
2460 : END DO
2461 8 : CALL cp_cfm_trace(tmp_conj_cfm, ref_state%overlap_cfm_kp(ikp), ztrace)
2462 18 : norm_ref = norm_ref + weight*REAL(ztrace, KIND=dp)
2463 : END DO
2464 2 : IF (PRESENT(para_env_inter_kp)) THEN
2465 2 : IF (ASSOCIATED(para_env_inter_kp)) THEN
2466 2 : CALL para_env_inter_kp%sum(error)
2467 2 : CALL para_env_inter_kp%sum(norm_ref)
2468 : END IF
2469 : END IF
2470 2 : IF (io_unit > 0) THEN
2471 1 : WRITE (UNIT=io_unit, FMT="(/,T2,A,F20.10)") "GEXT overlap fitting error:", &
2472 2 : SQRT(error/MAX(norm_ref, TINY(1.0_dp)))
2473 : END IF
2474 : END IF
2475 :
2476 2 : CALL cp_cfm_release(target_diff_cfm)
2477 2 : CALL cp_cfm_release(tmp_i_cfm)
2478 2 : CALL cp_cfm_release(tmp_j_cfm)
2479 2 : CALL cp_cfm_release(tmp_conj_cfm)
2480 2 : DEALLOCATE (A, b)
2481 2 : RETURN
2482 : END IF
2483 :
2484 696 : ALLOCATE (A(nvec - 1, nvec - 1), b(nvec - 1))
2485 :
2486 : ! get the reference for the difference fitting
2487 116 : ref_state => wfi_get_snapshot(wf_history, wf_index=1)
2488 :
2489 : ! assemble the target difference
2490 116 : CALL dbcsr_copy(target_diff, current_overlap)
2491 116 : CALL dbcsr_add(target_diff, ref_state%overlap, 1.0_dp, -1.0_dp)
2492 :
2493 : ! allocate tmp_k
2494 116 : CALL dbcsr_copy(tmp_k, current_overlap)
2495 :
2496 : ! assemble the matrix A and the RHS b
2497 348 : DO i = 2, nvec
2498 232 : state => wfi_get_snapshot(wf_history, wf_index=i)
2499 232 : CALL dbcsr_copy(tmp_i, state%overlap)
2500 232 : CALL dbcsr_add(tmp_i, ref_state%overlap, 1.0_dp, -1.0_dp)
2501 232 : CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_i, target_diff, 0.0_dp, tmp_k)
2502 232 : CALL dbcsr_trace(tmp_k, b(i - 1))
2503 :
2504 724 : DO j = 2, i
2505 376 : state => wfi_get_snapshot(wf_history, wf_index=j)
2506 376 : CALL dbcsr_copy(tmp_j, state%overlap)
2507 376 : CALL dbcsr_add(tmp_j, ref_state%overlap, 1.0_dp, -1.0_dp)
2508 376 : CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_i, tmp_j, 0.0_dp, tmp_k)
2509 376 : CALL dbcsr_trace(tmp_k, A(j - 1, i - 1))
2510 608 : A(i - 1, j - 1) = A(j - 1, i - 1)
2511 : END DO
2512 : END DO
2513 :
2514 : ! add the Tikhonov regularization
2515 348 : DO i = 1, nvec - 1
2516 348 : A(i, i) = A(i, i) + eps**2
2517 : END DO
2518 :
2519 : ! solve the linear system
2520 116 : CALL dposv('u', nvec - 1, 1, A, nvec - 1, b, nvec - 1, info)
2521 116 : IF (info /= 0) THEN
2522 0 : CPABORT("DPOSV failed.")
2523 : END IF
2524 :
2525 : ! set the coefficient for the reference snapshot
2526 348 : coeffs(1) = 1.0_dp - SUM(b)
2527 348 : 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 116 : IF (print_level > low_print_level) THEN
2532 20 : CALL dbcsr_copy(tmp_i, current_overlap)
2533 96 : DO i = 1, nvec
2534 76 : state => wfi_get_snapshot(wf_history, wf_index=i)
2535 96 : CALL dbcsr_add(tmp_i, state%overlap, 1.0_dp, -coeffs(i))
2536 : END DO
2537 20 : error = dbcsr_frobenius_norm(tmp_i)/dbcsr_frobenius_norm(state%overlap)
2538 20 : IF (io_unit > 0) THEN
2539 10 : 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 116 : CALL dbcsr_release(tmp_i)
2545 116 : CALL dbcsr_release(tmp_j)
2546 116 : CALL dbcsr_release(tmp_k)
2547 116 : CALL dbcsr_release(target_diff)
2548 116 : DEALLOCATE (A, b)
2549 :
2550 194 : 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 28 : SUBROUTINE tr_fitting(wf_history, current_overlap, coeffs, nvec, eps, io_unit, print_level, &
2573 28 : 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 28 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: b
2590 28 : 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 28 : IF (nvec <= 0) THEN
2597 0 : CPABORT("Not enough vectors to do the fitting")
2598 28 : ELSE IF (nvec == 1) THEN
2599 6 : coeffs(1) = 1.0_dp
2600 22 : RETURN
2601 : END IF
2602 :
2603 22 : IF (MOD(nvec, 2) == 0) THEN
2604 10 : ntr = nvec/2
2605 : ELSE
2606 12 : ntr = (nvec - 1)/2
2607 : END IF
2608 :
2609 22 : IF (PRESENT(current_overlap_kp)) THEN
2610 12 : ALLOCATE (A(ntr, ntr), b(ntr))
2611 2 : A = 0.0_dp
2612 2 : b = 0.0_dp
2613 :
2614 2 : ref_state => wfi_get_snapshot(wf_history, wf_index=nvec)
2615 2 : CALL cp_cfm_create(target_overlap_cfm, current_overlap_kp(1)%matrix_struct)
2616 2 : CALL cp_cfm_create(tmp_i_cfm, current_overlap_kp(1)%matrix_struct)
2617 2 : CALL cp_cfm_create(tmp_j_cfm, current_overlap_kp(1)%matrix_struct)
2618 2 : CALL cp_cfm_create(tmp_conj_cfm, current_overlap_kp(1)%matrix_struct)
2619 :
2620 10 : DO ikp = 1, SIZE(current_overlap_kp)
2621 8 : weight = 1.0_dp
2622 8 : IF (PRESENT(kpoint_weights)) weight = kpoint_weights(ikp)
2623 :
2624 8 : CALL cp_cfm_to_cfm(current_overlap_kp(ikp), target_overlap_cfm)
2625 8 : CALL cp_cfm_scale_and_add(z_one, target_overlap_cfm, z_one, ref_state%overlap_cfm_kp(ikp))
2626 18 : DO i = 1, ntr
2627 8 : state => wfi_get_snapshot(wf_history, wf_index=i)
2628 8 : CALL cp_cfm_to_cfm(state%overlap_cfm_kp(ikp), tmp_i_cfm)
2629 8 : state => wfi_get_snapshot(wf_history, wf_index=nvec - i)
2630 8 : CALL cp_cfm_scale_and_add(z_one, tmp_i_cfm, z_one, state%overlap_cfm_kp(ikp))
2631 :
2632 8 : CALL cp_cfm_to_cfm(tmp_i_cfm, tmp_conj_cfm)
2633 264 : DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
2634 4360 : DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
2635 : tmp_conj_cfm%local_data(irow_local, icol_local) = &
2636 4352 : CONJG(tmp_conj_cfm%local_data(irow_local, icol_local))
2637 : END DO
2638 : END DO
2639 8 : CALL cp_cfm_trace(tmp_conj_cfm, target_overlap_cfm, ztrace)
2640 8 : b(i) = b(i) + weight*REAL(ztrace, KIND=dp)
2641 24 : DO j = 1, i
2642 8 : state => wfi_get_snapshot(wf_history, wf_index=j)
2643 8 : CALL cp_cfm_to_cfm(state%overlap_cfm_kp(ikp), tmp_j_cfm)
2644 8 : state => wfi_get_snapshot(wf_history, wf_index=nvec - j)
2645 8 : CALL cp_cfm_scale_and_add(z_one, tmp_j_cfm, z_one, state%overlap_cfm_kp(ikp))
2646 8 : CALL cp_cfm_trace(tmp_conj_cfm, tmp_j_cfm, ztrace)
2647 16 : A(j, i) = A(j, i) + weight*REAL(ztrace, KIND=dp)
2648 : END DO
2649 : END DO
2650 : END DO
2651 :
2652 4 : DO i = 1, ntr
2653 6 : DO j = 1, i
2654 4 : A(i, j) = A(j, i)
2655 : END DO
2656 : END DO
2657 :
2658 2 : IF (PRESENT(para_env_inter_kp)) THEN
2659 2 : IF (ASSOCIATED(para_env_inter_kp)) THEN
2660 2 : CALL para_env_inter_kp%sum(A)
2661 2 : CALL para_env_inter_kp%sum(b)
2662 : END IF
2663 : END IF
2664 :
2665 4 : DO i = 1, ntr
2666 4 : A(i, i) = A(i, i) + eps**2
2667 : END DO
2668 :
2669 2 : CALL dposv('u', ntr, 1, A, ntr, b, ntr, info)
2670 2 : IF (info /= 0) THEN
2671 0 : CPABORT("DPOSV failed.")
2672 : END IF
2673 :
2674 6 : coeffs = 0.0_dp
2675 2 : coeffs(nvec) = -1.0_dp
2676 4 : DO i = 1, ntr
2677 2 : coeffs(i) = coeffs(i) + b(i)
2678 4 : coeffs(nvec - i) = coeffs(nvec - i) + b(i)
2679 : END DO
2680 :
2681 2 : IF (print_level > low_print_level) THEN
2682 2 : error = 0.0_dp
2683 2 : norm_ref = 0.0_dp
2684 10 : DO ikp = 1, SIZE(current_overlap_kp)
2685 8 : weight = 1.0_dp
2686 8 : IF (PRESENT(kpoint_weights)) weight = kpoint_weights(ikp)
2687 8 : CALL cp_cfm_to_cfm(current_overlap_kp(ikp), tmp_i_cfm)
2688 24 : DO i = 1, nvec
2689 16 : 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 24 : state%overlap_cfm_kp(ikp))
2692 : END DO
2693 8 : CALL cp_cfm_to_cfm(tmp_i_cfm, tmp_conj_cfm)
2694 264 : DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
2695 4360 : DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
2696 : tmp_conj_cfm%local_data(irow_local, icol_local) = &
2697 4352 : CONJG(tmp_conj_cfm%local_data(irow_local, icol_local))
2698 : END DO
2699 : END DO
2700 8 : CALL cp_cfm_trace(tmp_conj_cfm, tmp_i_cfm, ztrace)
2701 8 : error = error + weight*REAL(ztrace, KIND=dp)
2702 8 : CALL cp_cfm_to_cfm(ref_state%overlap_cfm_kp(ikp), tmp_conj_cfm)
2703 264 : DO icol_local = 1, SIZE(tmp_conj_cfm%local_data, 2)
2704 4360 : DO irow_local = 1, SIZE(tmp_conj_cfm%local_data, 1)
2705 : tmp_conj_cfm%local_data(irow_local, icol_local) = &
2706 4352 : CONJG(tmp_conj_cfm%local_data(irow_local, icol_local))
2707 : END DO
2708 : END DO
2709 8 : CALL cp_cfm_trace(tmp_conj_cfm, ref_state%overlap_cfm_kp(ikp), ztrace)
2710 18 : norm_ref = norm_ref + weight*REAL(ztrace, KIND=dp)
2711 : END DO
2712 2 : IF (PRESENT(para_env_inter_kp)) THEN
2713 2 : IF (ASSOCIATED(para_env_inter_kp)) THEN
2714 2 : CALL para_env_inter_kp%sum(error)
2715 2 : CALL para_env_inter_kp%sum(norm_ref)
2716 : END IF
2717 : END IF
2718 2 : IF (io_unit > 0) THEN
2719 1 : WRITE (UNIT=io_unit, FMT="(/,T2,A,F20.10)") "GEXT overlap fitting error:", &
2720 2 : SQRT(error/MAX(norm_ref, TINY(1.0_dp)))
2721 : END IF
2722 : END IF
2723 :
2724 2 : CALL cp_cfm_release(target_overlap_cfm)
2725 2 : CALL cp_cfm_release(tmp_i_cfm)
2726 2 : CALL cp_cfm_release(tmp_j_cfm)
2727 2 : CALL cp_cfm_release(tmp_conj_cfm)
2728 2 : DEALLOCATE (A, b)
2729 2 : RETURN
2730 : END IF
2731 :
2732 120 : ALLOCATE (A(ntr, ntr), b(ntr))
2733 :
2734 : ! get the reference for the difference fitting
2735 20 : ref_state => wfi_get_snapshot(wf_history, wf_index=nvec)
2736 :
2737 : ! assemble the target sum
2738 20 : CALL dbcsr_copy(target_overlap, current_overlap)
2739 20 : CALL dbcsr_add(target_overlap, ref_state%overlap, 1.0_dp, 1.0_dp)
2740 :
2741 : ! allocate tmp_k
2742 20 : CALL dbcsr_copy(tmp_k, current_overlap)
2743 :
2744 : ! assemble the matrix A and the RHS b
2745 52 : DO i = 1, ntr
2746 32 : state => wfi_get_snapshot(wf_history, wf_index=i)
2747 32 : CALL dbcsr_copy(tmp_i, state%overlap)
2748 32 : state => wfi_get_snapshot(wf_history, wf_index=nvec - i)
2749 32 : 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 32 : 0.0_dp, tmp_k)
2753 32 : CALL dbcsr_trace(tmp_k, b(i))
2754 96 : DO j = 1, i
2755 44 : state => wfi_get_snapshot(wf_history, wf_index=j)
2756 44 : CALL dbcsr_copy(tmp_j, state%overlap)
2757 44 : state => wfi_get_snapshot(wf_history, wf_index=nvec - j)
2758 44 : CALL dbcsr_add(tmp_j, state%overlap, 1.0_dp, 1.0_dp)
2759 44 : CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_i, tmp_j, 0.0_dp, tmp_k)
2760 44 : CALL dbcsr_trace(tmp_k, A(j, i))
2761 76 : A(i, j) = A(j, i)
2762 : END DO
2763 : END DO
2764 :
2765 : ! add the Tikhonov regularization
2766 52 : DO i = 1, ntr
2767 52 : A(i, i) = A(i, i) + eps**2
2768 : END DO
2769 :
2770 : ! solve the linear system
2771 20 : CALL dposv('u', ntr, 1, A, ntr, b, ntr, info)
2772 20 : IF (info /= 0) THEN
2773 0 : CPABORT("DPOSV failed.")
2774 : END IF
2775 :
2776 : ! reorder the coefficients
2777 96 : coeffs = 0.0_dp
2778 20 : coeffs(nvec) = -1.0_dp
2779 52 : DO i = 1, ntr
2780 32 : coeffs(i) = coeffs(i) + b(i)
2781 52 : 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 20 : IF (print_level > low_print_level) THEN
2787 20 : CALL dbcsr_copy(tmp_i, current_overlap)
2788 96 : DO i = 1, nvec
2789 76 : state => wfi_get_snapshot(wf_history, wf_index=i)
2790 96 : CALL dbcsr_add(tmp_i, state%overlap, 1.0_dp, -coeffs(i))
2791 : END DO
2792 20 : error = dbcsr_frobenius_norm(tmp_i)/dbcsr_frobenius_norm(state%overlap)
2793 20 : IF (io_unit > 0) THEN
2794 10 : 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 20 : CALL dbcsr_release(tmp_i)
2800 20 : CALL dbcsr_release(tmp_j)
2801 20 : CALL dbcsr_release(tmp_k)
2802 20 : CALL dbcsr_release(target_overlap)
2803 20 : DEALLOCATE (A, b)
2804 :
2805 44 : END SUBROUTINE tr_fitting
2806 :
2807 : END MODULE qs_wf_history_methods
|
