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