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