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