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 cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale,&
27 : cp_cfm_gemm,&
28 : cp_cfm_scale_and_add,&
29 : cp_cfm_scale_and_add_fm,&
30 : cp_cfm_triangular_multiply
31 : USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose
32 : USE cp_cfm_diag, ONLY: cp_cfm_heevd
33 : USE cp_cfm_types, ONLY: cp_cfm_create,&
34 : cp_cfm_release,&
35 : cp_cfm_set_all,&
36 : cp_cfm_to_cfm,&
37 : cp_cfm_to_fm,&
38 : cp_cfm_type,&
39 : 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_p_type, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
44 : dbcsr_type_symmetric
45 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
46 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
47 : cp_dbcsr_sm_fm_multiply,&
48 : dbcsr_allocate_matrix_set,&
49 : dbcsr_deallocate_matrix_set
50 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale,&
51 : cp_fm_scale_and_add
52 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
53 : cp_fm_pool_type,&
54 : fm_pool_create_fm,&
55 : fm_pool_give_back_fm,&
56 : fm_pools_create_fm_vect,&
57 : fm_pools_give_back_fm_vect
58 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
59 : cp_fm_struct_release,&
60 : cp_fm_struct_type
61 : USE cp_fm_types, ONLY: &
62 : copy_info_type, cp_fm_cleanup_copy_general, cp_fm_create, cp_fm_finish_copy_general, &
63 : cp_fm_get_info, cp_fm_release, cp_fm_start_copy_general, cp_fm_to_fm, cp_fm_type
64 : USE cp_log_handling, ONLY: cp_get_default_logger,&
65 : cp_logger_type,&
66 : cp_to_string
67 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
68 : cp_print_key_unit_nr,&
69 : low_print_level
70 : USE input_constants, ONLY: &
71 : wfi_aspc_nr, wfi_frozen_method_nr, wfi_linear_p_method_nr, wfi_linear_ps_method_nr, &
72 : wfi_linear_wf_method_nr, wfi_ps_method_nr, wfi_use_guess_method_nr, &
73 : wfi_use_prev_p_method_nr, wfi_use_prev_rho_r_method_nr, wfi_use_prev_wf_method_nr
74 : USE kinds, ONLY: dp
75 : USE kpoint_methods, ONLY: kpoint_density_matrices,&
76 : kpoint_density_transform,&
77 : kpoint_set_mo_occupation,&
78 : rskp_transform
79 : USE kpoint_types, ONLY: get_kpoint_info,&
80 : kpoint_env_type,&
81 : kpoint_type
82 : USE mathconstants, ONLY: gaussi,&
83 : z_one,&
84 : z_zero
85 : USE mathlib, ONLY: binomial
86 : USE message_passing, ONLY: mp_para_env_type
87 : USE parallel_gemm_api, ONLY: parallel_gemm
88 : USE pw_env_types, ONLY: pw_env_get,&
89 : pw_env_type
90 : USE pw_methods, ONLY: pw_copy
91 : USE pw_pool_types, ONLY: pw_pool_type
92 : USE pw_types, ONLY: pw_c1d_gs_type,&
93 : pw_r3d_rs_type
94 : USE qs_density_matrices, ONLY: calculate_density_matrix
95 : USE qs_environment_types, ONLY: get_qs_env,&
96 : qs_environment_type,&
97 : set_qs_env
98 : USE qs_ks_types, ONLY: qs_ks_did_change
99 : USE qs_matrix_pools, ONLY: mpools_get,&
100 : qs_matrix_pools_type
101 : USE qs_mo_methods, ONLY: make_basis_cholesky,&
102 : make_basis_lowdin,&
103 : make_basis_simple,&
104 : make_basis_sm
105 : USE qs_mo_types, ONLY: get_mo_set,&
106 : mo_set_type
107 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
108 : USE qs_rho_methods, ONLY: qs_rho_update_rho
109 : USE qs_rho_types, ONLY: qs_rho_get,&
110 : qs_rho_set,&
111 : qs_rho_type
112 : USE qs_scf_types, ONLY: ot_method_nr,&
113 : qs_scf_env_type
114 : USE qs_wf_history_types, ONLY: qs_wf_history_type,&
115 : qs_wf_snapshot_type,&
116 : wfi_get_snapshot,&
117 : wfi_release
118 : USE scf_control_types, ONLY: scf_control_type
119 : #include "./base/base_uses.f90"
120 :
121 : IMPLICIT NONE
122 : PRIVATE
123 :
124 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
125 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_wf_history_methods'
126 :
127 : PUBLIC :: wfi_create, wfi_update, wfi_create_for_kp, &
128 : wfi_extrapolate, wfi_get_method_label, &
129 : reorthogonalize_vectors, wfi_purge_history
130 :
131 : CONTAINS
132 :
133 : ! **************************************************************************************************
134 : !> \brief allocates and initialize a wavefunction snapshot
135 : !> \param snapshot the snapshot to create
136 : !> \par History
137 : !> 02.2003 created [fawzi]
138 : !> 02.2005 added wf_mol [MI]
139 : !> \author fawzi
140 : ! **************************************************************************************************
141 11334 : SUBROUTINE wfs_create(snapshot)
142 : TYPE(qs_wf_snapshot_type), INTENT(OUT) :: snapshot
143 :
144 : NULLIFY (snapshot%wf, snapshot%rho_r, &
145 : snapshot%rho_g, snapshot%rho_ao, snapshot%rho_ao_kp, &
146 : snapshot%overlap, snapshot%wf_kp, snapshot%overlap_cfm_kp, &
147 : snapshot%rho_frozen)
148 11334 : snapshot%dt = 1.0_dp
149 11334 : END SUBROUTINE wfs_create
150 :
151 : ! **************************************************************************************************
152 : !> \brief updates the given snapshot
153 : !> \param snapshot the snapshot to be updated
154 : !> \param wf_history the history
155 : !> \param qs_env the qs_env that should be snapshotted
156 : !> \param dt the time of the snapshot (wrt. to the previous snapshot)
157 : !> \par History
158 : !> 02.2003 created [fawzi]
159 : !> 02.2005 added kg_fm_mol_set for KG_GPW [MI]
160 : !> \author fawzi
161 : ! **************************************************************************************************
162 20386 : SUBROUTINE wfs_update(snapshot, wf_history, qs_env, dt)
163 : TYPE(qs_wf_snapshot_type), POINTER :: snapshot
164 : TYPE(qs_wf_history_type), POINTER :: wf_history
165 : TYPE(qs_environment_type), POINTER :: qs_env
166 : REAL(KIND=dp), INTENT(in), OPTIONAL :: dt
167 :
168 : CHARACTER(len=*), PARAMETER :: routineN = 'wfs_update'
169 :
170 : INTEGER :: handle, ic, igroup, ik, ikp, img, &
171 : indx_ft, ispin, kplocal, nc, nimg, &
172 : nkp_all, nkp_grps, nspin_kp, nspins
173 : INTEGER, DIMENSION(2) :: kp_range
174 20386 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
175 20386 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
176 : LOGICAL :: my_kpgrp
177 20386 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
178 20386 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info_ft
179 20386 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools, ao_mo_pools
180 : TYPE(cp_fm_struct_type), POINTER :: ao_ao_struct_ft
181 : TYPE(cp_fm_type) :: fmdummy_ft, fmlocal_ft
182 : TYPE(cp_fm_type), POINTER :: mo_coeff
183 20386 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
184 20386 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, rho_ao_kp
185 : TYPE(dbcsr_type), POINTER :: cmat_ft, rmat_ft, tmpmat_ft
186 : TYPE(dft_control_type), POINTER :: dft_control
187 : TYPE(kpoint_env_type), POINTER :: kp
188 : TYPE(kpoint_type), POINTER :: kpoints
189 20386 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
190 : TYPE(mp_para_env_type), POINTER :: para_env_ft
191 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
192 20386 : POINTER :: sab_nl
193 20386 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
194 : TYPE(pw_env_type), POINTER :: pw_env
195 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
196 20386 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
197 : TYPE(qs_matrix_pools_type), POINTER :: mpools_kp
198 : TYPE(qs_rho_type), POINTER :: rho
199 : TYPE(qs_scf_env_type), POINTER :: scf_env
200 :
201 20386 : CALL timeset(routineN, handle)
202 :
203 20386 : NULLIFY (pw_env, auxbas_pw_pool, ao_mo_pools, ao_ao_fm_pools, dft_control, mos, mo_coeff, &
204 20386 : rho, rho_r, rho_g, rho_ao, matrix_s, matrix_s_kp, kpoints, kp, &
205 20386 : kp_dist, cell_to_index, xkp, sab_nl, scf_env, mpools_kp, para_env_ft, &
206 20386 : rmat_ft, cmat_ft, tmpmat_ft, ao_ao_struct_ft)
207 : CALL get_qs_env(qs_env, pw_env=pw_env, &
208 20386 : dft_control=dft_control, rho=rho)
209 20386 : CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_pools)
210 20386 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
211 :
212 20386 : CPASSERT(ASSOCIATED(wf_history))
213 20386 : CPASSERT(ASSOCIATED(dft_control))
214 20386 : IF (.NOT. ASSOCIATED(snapshot)) THEN
215 11334 : ALLOCATE (snapshot)
216 11334 : CALL wfs_create(snapshot)
217 : END IF
218 20386 : CPASSERT(wf_history%ref_count > 0)
219 :
220 20386 : nspins = dft_control%nspins
221 20386 : snapshot%dt = 1.0_dp
222 20386 : IF (PRESENT(dt)) snapshot%dt = dt
223 20386 : IF (wf_history%store_wf) THEN
224 19340 : CALL get_qs_env(qs_env, mos=mos)
225 19340 : IF (.NOT. ASSOCIATED(snapshot%wf)) THEN
226 : CALL fm_pools_create_fm_vect(ao_mo_pools, snapshot%wf, &
227 11056 : name="ws_snap-ws")
228 11056 : CPASSERT(nspins == SIZE(snapshot%wf))
229 : END IF
230 40804 : DO ispin = 1, nspins
231 21464 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
232 40804 : CALL cp_fm_to_fm(mo_coeff, snapshot%wf(ispin))
233 : END DO
234 : ELSE
235 1046 : CALL fm_pools_give_back_fm_vect(ao_mo_pools, snapshot%wf)
236 : END IF
237 :
238 20386 : IF (wf_history%store_rho_r) THEN
239 0 : CALL qs_rho_get(rho, rho_r=rho_r)
240 0 : CPASSERT(ASSOCIATED(rho_r))
241 0 : IF (.NOT. ASSOCIATED(snapshot%rho_r)) THEN
242 0 : ALLOCATE (snapshot%rho_r(nspins))
243 0 : DO ispin = 1, nspins
244 0 : CALL auxbas_pw_pool%create_pw(snapshot%rho_r(ispin))
245 : END DO
246 : END IF
247 0 : DO ispin = 1, nspins
248 0 : CALL pw_copy(rho_r(ispin), snapshot%rho_r(ispin))
249 : END DO
250 20386 : ELSE IF (ASSOCIATED(snapshot%rho_r)) THEN
251 0 : DO ispin = 1, SIZE(snapshot%rho_r)
252 0 : CALL auxbas_pw_pool%give_back_pw(snapshot%rho_r(ispin))
253 : END DO
254 0 : DEALLOCATE (snapshot%rho_r)
255 : END IF
256 :
257 20386 : IF (wf_history%store_rho_g) THEN
258 0 : CALL qs_rho_get(rho, rho_g=rho_g)
259 0 : CPASSERT(ASSOCIATED(rho_g))
260 0 : IF (.NOT. ASSOCIATED(snapshot%rho_g)) THEN
261 0 : ALLOCATE (snapshot%rho_g(nspins))
262 0 : DO ispin = 1, nspins
263 0 : CALL auxbas_pw_pool%create_pw(snapshot%rho_g(ispin))
264 : END DO
265 : END IF
266 0 : DO ispin = 1, nspins
267 0 : CALL pw_copy(rho_g(ispin), snapshot%rho_g(ispin))
268 : END DO
269 20386 : ELSE IF (ASSOCIATED(snapshot%rho_g)) THEN
270 0 : DO ispin = 1, SIZE(snapshot%rho_g)
271 0 : CALL auxbas_pw_pool%give_back_pw(snapshot%rho_g(ispin))
272 : END DO
273 0 : DEALLOCATE (snapshot%rho_g)
274 : END IF
275 :
276 20386 : IF (ASSOCIATED(snapshot%rho_ao)) THEN ! the sparsity might be different
277 : ! (future struct:check)
278 270 : CALL dbcsr_deallocate_matrix_set(snapshot%rho_ao)
279 : END IF
280 20386 : IF (wf_history%store_rho_ao) THEN
281 338 : CALL qs_rho_get(rho, rho_ao=rho_ao)
282 338 : CPASSERT(ASSOCIATED(rho_ao))
283 :
284 338 : CALL dbcsr_allocate_matrix_set(snapshot%rho_ao, nspins)
285 836 : DO ispin = 1, nspins
286 498 : ALLOCATE (snapshot%rho_ao(ispin)%matrix)
287 836 : CALL dbcsr_copy(snapshot%rho_ao(ispin)%matrix, rho_ao(ispin)%matrix)
288 : END DO
289 : END IF
290 :
291 20386 : IF (ASSOCIATED(snapshot%rho_ao_kp)) THEN ! the sparsity might be different
292 : ! (future struct:check)
293 220 : CALL dbcsr_deallocate_matrix_set(snapshot%rho_ao_kp)
294 : END IF
295 20386 : IF (wf_history%store_rho_ao_kp) THEN
296 232 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
297 232 : CPASSERT(ASSOCIATED(rho_ao_kp))
298 :
299 232 : nimg = dft_control%nimages
300 232 : CALL dbcsr_allocate_matrix_set(snapshot%rho_ao_kp, nspins, nimg)
301 554 : DO ispin = 1, nspins
302 34092 : DO img = 1, nimg
303 33538 : ALLOCATE (snapshot%rho_ao_kp(ispin, img)%matrix)
304 : CALL dbcsr_copy(snapshot%rho_ao_kp(ispin, img)%matrix, &
305 33860 : rho_ao_kp(ispin, img)%matrix)
306 : END DO
307 : END DO
308 : END IF
309 :
310 20386 : IF (ASSOCIATED(snapshot%overlap)) THEN ! the sparsity might be different
311 : ! (future struct:check)
312 5676 : CALL dbcsr_deallocate_matrix(snapshot%overlap)
313 : END IF
314 20386 : IF (wf_history%store_overlap) THEN
315 15698 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
316 15698 : CPASSERT(ASSOCIATED(matrix_s))
317 15698 : CPASSERT(ASSOCIATED(matrix_s(1)%matrix))
318 15698 : ALLOCATE (snapshot%overlap)
319 15698 : CALL dbcsr_copy(snapshot%overlap, matrix_s(1)%matrix)
320 : END IF
321 :
322 20386 : CALL get_qs_env(qs_env, kpoints=kpoints)
323 20386 : IF (ASSOCIATED(kpoints)) THEN
324 20386 : IF (ASSOCIATED(kpoints%kp_env)) THEN
325 : ! --- k-point WFN snapshot: store complex MO coefficients per local k-point ---
326 692 : IF (wf_history%store_wf_kp) THEN
327 460 : CALL get_kpoint_info(kpoints, kp_range=kp_range)
328 460 : kplocal = kp_range(2) - kp_range(1) + 1
329 460 : nspin_kp = SIZE(kpoints%kp_env(1)%kpoint_env%mos, 2)
330 460 : nc = SIZE(kpoints%kp_env(1)%kpoint_env%mos, 1) ! 2=complex, 1=real
331 :
332 460 : IF (ASSOCIATED(snapshot%wf_kp)) THEN
333 732 : DO ikp = 1, SIZE(snapshot%wf_kp, 1)
334 1664 : DO ic = 1, SIZE(snapshot%wf_kp, 2)
335 2330 : DO ispin = 1, SIZE(snapshot%wf_kp, 3)
336 1864 : CALL cp_fm_release(snapshot%wf_kp(ikp, ic, ispin))
337 : END DO
338 : END DO
339 : END DO
340 266 : DEALLOCATE (snapshot%wf_kp)
341 : END IF
342 :
343 7388 : ALLOCATE (snapshot%wf_kp(kplocal, nc, nspin_kp))
344 1950 : DO ikp = 1, kplocal
345 1490 : kp => kpoints%kp_env(ikp)%kpoint_env
346 3780 : DO ispin = 1, nspin_kp
347 6980 : DO ic = 1, nc
348 3660 : CALL get_mo_set(kp%mos(ic, ispin), mo_coeff=mo_coeff)
349 : CALL cp_fm_create(snapshot%wf_kp(ikp, ic, ispin), &
350 : mo_coeff%matrix_struct, &
351 3660 : name="wfkp_snap")
352 5490 : CALL cp_fm_to_fm(mo_coeff, snapshot%wf_kp(ikp, ic, ispin))
353 : END DO
354 : END DO
355 : END DO
356 : END IF
357 :
358 : ! --- k-point overlap snapshot: Fourier-transform S(R)→S(k) NOW and store as cfm ---
359 : ! This is critical: we MUST transform at snapshot time using the CURRENT neighbor
360 : ! list. Storing S(R) and re-transforming later would use a stale neighbor list,
361 : ! producing wrong S(k) when the neighbor list changes during MD.
362 692 : IF (wf_history%store_overlap_kp) THEN
363 460 : CALL get_qs_env(qs_env, matrix_s_kp=matrix_s_kp, scf_env=scf_env)
364 : CALL get_kpoint_info(kpoints, nkp=nkp_all, xkp=xkp, kp_range=kp_range, &
365 : nkp_groups=nkp_grps, kp_dist=kp_dist, &
366 460 : sab_nl=sab_nl, cell_to_index=cell_to_index)
367 460 : kplocal = kp_range(2) - kp_range(1) + 1
368 460 : para_env_ft => kpoints%blacs_env_all%para_env
369 :
370 : ! Allocate dbcsr work matrices for FT (same pattern as do_general_diag_kp)
371 460 : ALLOCATE (rmat_ft, cmat_ft, tmpmat_ft)
372 : CALL dbcsr_create(rmat_ft, template=matrix_s_kp(1, 1)%matrix, &
373 460 : matrix_type=dbcsr_type_symmetric)
374 : CALL dbcsr_create(cmat_ft, template=matrix_s_kp(1, 1)%matrix, &
375 460 : matrix_type=dbcsr_type_antisymmetric)
376 : CALL dbcsr_create(tmpmat_ft, template=matrix_s_kp(1, 1)%matrix, &
377 460 : matrix_type=dbcsr_type_no_symmetry)
378 460 : CALL cp_dbcsr_alloc_block_from_nbl(rmat_ft, sab_nl)
379 460 : CALL cp_dbcsr_alloc_block_from_nbl(cmat_ft, sab_nl)
380 :
381 : ! Get kp-subgroup FM from pool
382 460 : CALL get_kpoint_info(kpoints, mpools=mpools_kp)
383 460 : CALL mpools_get(mpools_kp, ao_ao_fm_pools=ao_ao_fm_pools)
384 460 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal_ft)
385 :
386 : ! Release old snapshot if present
387 460 : IF (ASSOCIATED(snapshot%overlap_cfm_kp)) THEN
388 732 : DO ikp = 1, SIZE(snapshot%overlap_cfm_kp)
389 732 : CALL cp_cfm_release(snapshot%overlap_cfm_kp(ikp))
390 : END DO
391 266 : DEALLOCATE (snapshot%overlap_cfm_kp)
392 : END IF
393 2870 : ALLOCATE (snapshot%overlap_cfm_kp(kplocal))
394 :
395 460 : CALL cp_fm_get_info(fmlocal_ft, matrix_struct=ao_ao_struct_ft)
396 :
397 : ! Communication info array
398 11028 : ALLOCATE (info_ft(kplocal*nkp_grps, 2))
399 :
400 : ! Phase A: Start async FT + redistribute for each k-point
401 460 : indx_ft = 0
402 1950 : DO ikp = 1, kplocal
403 4474 : DO igroup = 1, nkp_grps
404 2524 : ik = kp_dist(1, igroup) + ikp - 1
405 2524 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
406 2524 : indx_ft = indx_ft + 1
407 :
408 2524 : CALL dbcsr_set(rmat_ft, 0.0_dp)
409 2524 : CALL dbcsr_set(cmat_ft, 0.0_dp)
410 : CALL rskp_transform(rmatrix=rmat_ft, cmatrix=cmat_ft, rsmat=matrix_s_kp, &
411 : ispin=1, xkp=xkp(1:3, ik), &
412 2524 : cell_to_index=cell_to_index, sab_nl=sab_nl)
413 2524 : CALL dbcsr_desymmetrize(rmat_ft, tmpmat_ft)
414 2524 : CALL copy_dbcsr_to_fm(tmpmat_ft, scf_env%scf_work1(1))
415 2524 : CALL dbcsr_desymmetrize(cmat_ft, tmpmat_ft)
416 2524 : CALL copy_dbcsr_to_fm(tmpmat_ft, scf_env%scf_work1(2))
417 :
418 4014 : IF (my_kpgrp) THEN
419 : CALL cp_fm_start_copy_general(scf_env%scf_work1(1), fmlocal_ft, &
420 1490 : para_env_ft, info_ft(indx_ft, 1))
421 : CALL cp_fm_start_copy_general(scf_env%scf_work1(2), fmlocal_ft, &
422 1490 : para_env_ft, info_ft(indx_ft, 2))
423 : ELSE
424 : CALL cp_fm_start_copy_general(scf_env%scf_work1(1), fmdummy_ft, &
425 1034 : para_env_ft, info_ft(indx_ft, 1))
426 : CALL cp_fm_start_copy_general(scf_env%scf_work1(2), fmdummy_ft, &
427 1034 : para_env_ft, info_ft(indx_ft, 2))
428 : END IF
429 : END DO
430 : END DO
431 :
432 : ! Phase B: Finish communication and assemble S(k) as cfm
433 : indx_ft = 0
434 1950 : DO ikp = 1, kplocal
435 1490 : CALL cp_cfm_create(snapshot%overlap_cfm_kp(ikp), ao_ao_struct_ft)
436 1490 : CALL cp_cfm_set_all(snapshot%overlap_cfm_kp(ikp), z_zero)
437 4474 : DO igroup = 1, nkp_grps
438 2524 : ik = kp_dist(1, igroup) + ikp - 1
439 2524 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
440 1034 : indx_ft = indx_ft + 1
441 1490 : IF (my_kpgrp) THEN
442 1490 : CALL cp_fm_finish_copy_general(fmlocal_ft, info_ft(indx_ft, 1))
443 : CALL cp_cfm_scale_and_add_fm(z_zero, snapshot%overlap_cfm_kp(ikp), &
444 1490 : z_one, fmlocal_ft)
445 1490 : CALL cp_fm_finish_copy_general(fmlocal_ft, info_ft(indx_ft, 2))
446 : CALL cp_cfm_scale_and_add_fm(z_one, snapshot%overlap_cfm_kp(ikp), &
447 1490 : gaussi, fmlocal_ft)
448 : END IF
449 : END DO
450 : END DO
451 :
452 : ! Cleanup
453 2984 : DO indx_ft = 1, kplocal*nkp_grps
454 2524 : CALL cp_fm_cleanup_copy_general(info_ft(indx_ft, 1))
455 2984 : CALL cp_fm_cleanup_copy_general(info_ft(indx_ft, 2))
456 : END DO
457 5968 : DEALLOCATE (info_ft)
458 460 : CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal_ft)
459 460 : CALL dbcsr_deallocate_matrix(rmat_ft)
460 460 : CALL dbcsr_deallocate_matrix(cmat_ft)
461 920 : CALL dbcsr_deallocate_matrix(tmpmat_ft)
462 : END IF
463 : END IF
464 : END IF
465 :
466 : IF (wf_history%store_frozen_density) THEN
467 : ! do nothing
468 : ! CALL deallocate_matrix_set(snapshot%rho_frozen%rho_ao)
469 : END IF
470 :
471 20386 : CALL timestop(handle)
472 :
473 20386 : END SUBROUTINE wfs_update
474 :
475 : ! **************************************************************************************************
476 : !> \brief ...
477 : !> \param wf_history ...
478 : !> \param interpolation_method_nr the tag of the method used for
479 : !> the extrapolation of the initial density for the next md step
480 : !> (see qs_wf_history_types:wfi_*_method_nr)
481 : !> \param extrapolation_order ...
482 : !> \param has_unit_metric ...
483 : !> \par History
484 : !> 02.2003 created [fawzi]
485 : !> \author fawzi
486 : ! **************************************************************************************************
487 7710 : SUBROUTINE wfi_create(wf_history, interpolation_method_nr, extrapolation_order, &
488 : has_unit_metric)
489 : TYPE(qs_wf_history_type), POINTER :: wf_history
490 : INTEGER, INTENT(in) :: interpolation_method_nr, &
491 : extrapolation_order
492 : LOGICAL, INTENT(IN) :: has_unit_metric
493 :
494 : CHARACTER(len=*), PARAMETER :: routineN = 'wfi_create'
495 :
496 : INTEGER :: handle, i
497 :
498 7710 : CALL timeset(routineN, handle)
499 :
500 7710 : ALLOCATE (wf_history)
501 7710 : wf_history%ref_count = 1
502 7710 : wf_history%memory_depth = 0
503 7710 : wf_history%snapshot_count = 0
504 7710 : wf_history%last_state_index = 1
505 : wf_history%store_wf = .FALSE.
506 : wf_history%store_rho_r = .FALSE.
507 : wf_history%store_rho_g = .FALSE.
508 : wf_history%store_rho_ao = .FALSE.
509 : wf_history%store_rho_ao_kp = .FALSE.
510 : wf_history%store_overlap = .FALSE.
511 : wf_history%store_wf_kp = .FALSE.
512 : wf_history%store_overlap_kp = .FALSE.
513 : wf_history%store_frozen_density = .FALSE.
514 : NULLIFY (wf_history%past_states)
515 :
516 7710 : wf_history%interpolation_method_nr = interpolation_method_nr
517 :
518 : SELECT CASE (wf_history%interpolation_method_nr)
519 : CASE (wfi_use_guess_method_nr)
520 : wf_history%memory_depth = 0
521 : CASE (wfi_use_prev_wf_method_nr)
522 64 : wf_history%memory_depth = 0
523 : CASE (wfi_use_prev_p_method_nr)
524 64 : wf_history%memory_depth = 1
525 64 : wf_history%store_rho_ao = .TRUE.
526 : CASE (wfi_use_prev_rho_r_method_nr)
527 4 : wf_history%memory_depth = 1
528 4 : wf_history%store_rho_ao = .TRUE.
529 : CASE (wfi_linear_wf_method_nr)
530 4 : wf_history%memory_depth = 2
531 4 : wf_history%store_wf = .TRUE.
532 : CASE (wfi_linear_p_method_nr)
533 6 : wf_history%memory_depth = 2
534 6 : wf_history%store_rho_ao = .TRUE.
535 : CASE (wfi_linear_ps_method_nr)
536 6 : wf_history%memory_depth = 2
537 6 : wf_history%store_wf = .TRUE.
538 6 : IF (.NOT. has_unit_metric) wf_history%store_overlap = .TRUE.
539 : CASE (wfi_ps_method_nr)
540 343 : CALL cite_reference(VandeVondele2005a)
541 343 : wf_history%memory_depth = extrapolation_order + 1
542 343 : wf_history%store_wf = .TRUE.
543 343 : wf_history%store_wf_kp = .TRUE.
544 343 : IF (.NOT. has_unit_metric) THEN
545 339 : wf_history%store_overlap = .TRUE.
546 339 : wf_history%store_overlap_kp = .TRUE.
547 : END IF
548 : CASE (wfi_frozen_method_nr)
549 4 : wf_history%memory_depth = 1
550 4 : wf_history%store_frozen_density = .TRUE.
551 : CASE (wfi_aspc_nr)
552 7141 : wf_history%memory_depth = extrapolation_order + 2
553 7141 : wf_history%store_wf = .TRUE.
554 7141 : wf_history%store_wf_kp = .TRUE.
555 7141 : IF (.NOT. has_unit_metric) THEN
556 6159 : wf_history%store_overlap = .TRUE.
557 6159 : wf_history%store_overlap_kp = .TRUE.
558 : END IF
559 : CASE default
560 : CALL cp_abort(__LOCATION__, &
561 : "Unknown interpolation method: "// &
562 0 : TRIM(ADJUSTL(cp_to_string(interpolation_method_nr))))
563 7710 : wf_history%interpolation_method_nr = wfi_use_prev_rho_r_method_nr
564 : END SELECT
565 59853 : ALLOCATE (wf_history%past_states(wf_history%memory_depth))
566 :
567 44571 : DO i = 1, SIZE(wf_history%past_states)
568 44571 : NULLIFY (wf_history%past_states(i)%snapshot)
569 : END DO
570 :
571 7710 : CALL timestop(handle)
572 7710 : END SUBROUTINE wfi_create
573 :
574 : ! **************************************************************************************************
575 : !> \brief Adapts wf_history storage flags for k-point calculations.
576 : !> For ASPC, switches from Gamma WFN storage to k-point WFN storage.
577 : !> Other WFN-based methods remain blocked.
578 : !> \param wf_history ...
579 : !> \par History
580 : !> 06.2015 created [jhu]
581 : !> \author jhu
582 : ! **************************************************************************************************
583 194 : SUBROUTINE wfi_create_for_kp(wf_history)
584 : TYPE(qs_wf_history_type), POINTER :: wf_history
585 :
586 194 : CPASSERT(ASSOCIATED(wf_history))
587 194 : IF (wf_history%store_rho_ao) THEN
588 10 : wf_history%store_rho_ao_kp = .TRUE.
589 10 : wf_history%store_rho_ao = .FALSE.
590 : END IF
591 : ! KP-compatible PS and ASPC: switch from Gamma WFN storage to KP WFN storage
592 194 : IF (wf_history%store_wf_kp) THEN
593 100 : wf_history%store_wf = .FALSE.
594 100 : wf_history%store_overlap = .FALSE.
595 : ! store_wf_kp and store_overlap_kp remain TRUE
596 : ELSE
597 : ! Linear methods (except LINEAR_P) are still blocked
598 94 : IF (wf_history%store_wf .OR. wf_history%store_overlap) THEN
599 0 : CPABORT("Linear WFN-based extrapolation methods not implemented for k-points.")
600 : END IF
601 : END IF
602 194 : IF (wf_history%store_frozen_density) THEN
603 0 : CPABORT("Frozen density initialization method not possible for kpoints.")
604 : END IF
605 :
606 194 : END SUBROUTINE wfi_create_for_kp
607 :
608 : ! **************************************************************************************************
609 : !> \brief returns a string describing the interpolation method
610 : !> \param method_nr ...
611 : !> \return ...
612 : !> \par History
613 : !> 02.2003 created [fawzi]
614 : !> \author fawzi
615 : ! **************************************************************************************************
616 10759 : FUNCTION wfi_get_method_label(method_nr) RESULT(res)
617 : INTEGER, INTENT(in) :: method_nr
618 : CHARACTER(len=30) :: res
619 :
620 10759 : res = "unknown"
621 10997 : SELECT CASE (method_nr)
622 : CASE (wfi_use_prev_p_method_nr)
623 238 : res = "previous_p"
624 : CASE (wfi_use_prev_wf_method_nr)
625 213 : res = "previous_wf"
626 : CASE (wfi_use_prev_rho_r_method_nr)
627 4 : res = "previous_rho_r"
628 : CASE (wfi_use_guess_method_nr)
629 3426 : res = "initial_guess"
630 : CASE (wfi_linear_wf_method_nr)
631 2 : res = "mo linear"
632 : CASE (wfi_linear_p_method_nr)
633 3 : res = "P linear"
634 : CASE (wfi_linear_ps_method_nr)
635 6 : res = "PS linear"
636 : CASE (wfi_ps_method_nr)
637 188 : res = "PS Nth order"
638 : CASE (wfi_frozen_method_nr)
639 4 : res = "frozen density approximation"
640 : CASE (wfi_aspc_nr)
641 6675 : res = "ASPC"
642 : CASE default
643 : CALL cp_abort(__LOCATION__, &
644 : "Unknown interpolation method: "// &
645 10759 : TRIM(ADJUSTL(cp_to_string(method_nr))))
646 : END SELECT
647 10759 : END FUNCTION wfi_get_method_label
648 :
649 : ! **************************************************************************************************
650 : !> \brief calculates the new starting state for the scf for the next
651 : !> wf optimization
652 : !> \param wf_history the previous history needed to extrapolate
653 : !> \param qs_env the qs env with the latest result, and that will contain
654 : !> the new starting state
655 : !> \param dt the time at which to extrapolate (wrt. to the last snapshot)
656 : !> \param extrapolation_method_nr returns the extrapolation method used
657 : !> \param orthogonal_wf ...
658 : !> \par History
659 : !> 02.2003 created [fawzi]
660 : !> 11.2003 Joost VandeVondele : Implemented Nth order PS extrapolation
661 : !> \author fawzi
662 : ! **************************************************************************************************
663 21169 : SUBROUTINE wfi_extrapolate(wf_history, qs_env, dt, extrapolation_method_nr, &
664 : orthogonal_wf)
665 : TYPE(qs_wf_history_type), POINTER :: wf_history
666 : TYPE(qs_environment_type), POINTER :: qs_env
667 : REAL(KIND=dp), INTENT(IN) :: dt
668 : INTEGER, INTENT(OUT), OPTIONAL :: extrapolation_method_nr
669 : LOGICAL, INTENT(OUT), OPTIONAL :: orthogonal_wf
670 :
671 : CHARACTER(len=*), PARAMETER :: routineN = 'wfi_extrapolate'
672 :
673 : INTEGER :: actual_extrapolation_method_nr, handle, &
674 : i, img, io_unit, ispin, k, n, nmo, &
675 : nvec, print_level
676 : LOGICAL :: do_kpoints, my_orthogonal_wf, use_overlap
677 : REAL(KIND=dp) :: alpha, t0, t1, t2
678 21169 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools
679 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct, matrix_struct_new
680 : TYPE(cp_fm_type) :: csc, fm_tmp
681 : TYPE(cp_fm_type), POINTER :: mo_coeff
682 : TYPE(cp_logger_type), POINTER :: logger
683 21169 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao, rho_frozen_ao
684 21169 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
685 21169 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
686 : TYPE(qs_rho_type), POINTER :: rho
687 : TYPE(qs_wf_snapshot_type), POINTER :: t0_state, t1_state
688 :
689 21169 : NULLIFY (mos, ao_mo_fm_pools, t0_state, t1_state, mo_coeff, &
690 21169 : rho, rho_ao, rho_frozen_ao)
691 :
692 21169 : use_overlap = wf_history%store_overlap
693 :
694 21169 : CALL timeset(routineN, handle)
695 21169 : logger => cp_get_default_logger()
696 21169 : print_level = logger%iter_info%print_level
697 : io_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%SCF%PRINT%PROGRAM_RUN_INFO", &
698 21169 : extension=".scfLog")
699 :
700 21169 : CPASSERT(ASSOCIATED(wf_history))
701 21169 : CPASSERT(wf_history%ref_count > 0)
702 21169 : CPASSERT(ASSOCIATED(qs_env))
703 21169 : CALL get_qs_env(qs_env, mos=mos, rho=rho, do_kpoints=do_kpoints)
704 21169 : CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
705 : ! chooses the method for this extrapolation
706 21169 : IF (wf_history%snapshot_count < 1) THEN
707 : actual_extrapolation_method_nr = wfi_use_guess_method_nr
708 : ELSE
709 14616 : actual_extrapolation_method_nr = wf_history%interpolation_method_nr
710 : END IF
711 :
712 8 : SELECT CASE (actual_extrapolation_method_nr)
713 : CASE (wfi_linear_wf_method_nr)
714 8 : IF (wf_history%snapshot_count < 2) THEN
715 4 : actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
716 : END IF
717 : CASE (wfi_linear_p_method_nr)
718 12 : IF (wf_history%snapshot_count < 2) THEN
719 6 : actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
720 : END IF
721 : CASE (wfi_linear_ps_method_nr)
722 14616 : IF (wf_history%snapshot_count < 2) THEN
723 6 : actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
724 : END IF
725 : END SELECT
726 :
727 21169 : IF (PRESENT(extrapolation_method_nr)) &
728 21169 : extrapolation_method_nr = actual_extrapolation_method_nr
729 21169 : my_orthogonal_wf = .FALSE.
730 :
731 8 : SELECT CASE (actual_extrapolation_method_nr)
732 : CASE (wfi_frozen_method_nr)
733 8 : CPASSERT(.NOT. do_kpoints)
734 8 : t0_state => wfi_get_snapshot(wf_history, wf_index=1)
735 8 : CPASSERT(ASSOCIATED(t0_state%rho_frozen))
736 :
737 8 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
738 8 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
739 :
740 8 : CALL qs_rho_get(t0_state%rho_frozen, rho_ao=rho_frozen_ao)
741 8 : CALL qs_rho_get(rho, rho_ao=rho_ao)
742 16 : DO ispin = 1, SIZE(rho_frozen_ao)
743 : CALL dbcsr_copy(rho_ao(ispin)%matrix, &
744 : rho_frozen_ao(ispin)%matrix, &
745 16 : keep_sparsity=.TRUE.)
746 : END DO
747 : !FM updating rho_ao directly with t0_state%rho_ao would have the
748 : !FM wrong matrix structure
749 8 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
750 8 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
751 :
752 8 : my_orthogonal_wf = .FALSE.
753 : CASE (wfi_use_prev_rho_r_method_nr)
754 8 : t0_state => wfi_get_snapshot(wf_history, wf_index=1)
755 8 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
756 8 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
757 8 : IF (do_kpoints) THEN
758 0 : CPASSERT(ASSOCIATED(t0_state%rho_ao_kp))
759 0 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
760 0 : DO ispin = 1, SIZE(t0_state%rho_ao_kp, 1)
761 0 : DO img = 1, SIZE(t0_state%rho_ao_kp, 2)
762 0 : IF (img > SIZE(rho_ao_kp, 2)) THEN
763 0 : CPWARN("Change in cell neighborlist: might affect quality of initial guess")
764 : ELSE
765 : CALL dbcsr_copy(rho_ao_kp(ispin, img)%matrix, &
766 : t0_state%rho_ao_kp(ispin, img)%matrix, &
767 0 : keep_sparsity=.TRUE.)
768 : END IF
769 : END DO
770 : END DO
771 : ELSE
772 8 : CPASSERT(ASSOCIATED(t0_state%rho_ao))
773 8 : CALL qs_rho_get(rho, rho_ao=rho_ao)
774 16 : DO ispin = 1, SIZE(t0_state%rho_ao)
775 : CALL dbcsr_copy(rho_ao(ispin)%matrix, &
776 : t0_state%rho_ao(ispin)%matrix, &
777 16 : keep_sparsity=.TRUE.)
778 : END DO
779 : END IF
780 : ! Why is rho_g valid at this point ?
781 8 : CALL qs_rho_set(rho, rho_g_valid=.TRUE.)
782 :
783 : ! does nothing
784 : CASE (wfi_use_prev_wf_method_nr)
785 426 : my_orthogonal_wf = .TRUE.
786 426 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
787 426 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
788 :
789 426 : IF (do_kpoints) THEN
790 6 : CALL wfi_use_prev_wf_kp(qs_env, io_unit, print_level)
791 : ELSE
792 420 : CALL qs_rho_get(rho, rho_ao=rho_ao)
793 888 : DO ispin = 1, SIZE(mos)
794 468 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
795 468 : CALL reorthogonalize_vectors(qs_env, v_matrix=mo_coeff, n_col=nmo)
796 1356 : CALL calculate_density_matrix(mo_set=mos(ispin), density_matrix=rho_ao(ispin)%matrix)
797 : END DO
798 420 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
799 420 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
800 : END IF
801 :
802 : CASE (wfi_use_prev_p_method_nr)
803 476 : t0_state => wfi_get_snapshot(wf_history, wf_index=1)
804 476 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
805 476 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
806 476 : IF (do_kpoints) THEN
807 218 : CPASSERT(ASSOCIATED(t0_state%rho_ao_kp))
808 218 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
809 524 : DO ispin = 1, SIZE(t0_state%rho_ao_kp, 1)
810 31248 : DO img = 1, SIZE(t0_state%rho_ao_kp, 2)
811 31030 : IF (img > SIZE(rho_ao_kp, 2)) THEN
812 18 : CPWARN("Change in cell neighborlist: might affect quality of initial guess")
813 : ELSE
814 : CALL dbcsr_copy(rho_ao_kp(ispin, img)%matrix, &
815 : t0_state%rho_ao_kp(ispin, img)%matrix, &
816 30706 : keep_sparsity=.TRUE.)
817 : END IF
818 : END DO
819 : END DO
820 : ELSE
821 258 : CPASSERT(ASSOCIATED(t0_state%rho_ao))
822 258 : CALL qs_rho_get(rho, rho_ao=rho_ao)
823 646 : DO ispin = 1, SIZE(t0_state%rho_ao)
824 : CALL dbcsr_copy(rho_ao(ispin)%matrix, &
825 : t0_state%rho_ao(ispin)%matrix, &
826 646 : keep_sparsity=.TRUE.)
827 : END DO
828 : END IF
829 : !FM updating rho_ao directly with t0_state%rho_ao would have the
830 : !FM wrong matrix structure
831 476 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
832 476 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
833 :
834 : CASE (wfi_use_guess_method_nr)
835 : !FM more clean to do it here, but it
836 : !FM might need to read a file (restart) and thus globenv
837 : !FM I do not want globenv here, thus done by the caller
838 : !FM (btw. it also needs the eigensolver, and unless you relocate it
839 : !FM gives circular dependencies)
840 6795 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
841 6795 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
842 : CASE (wfi_linear_wf_method_nr)
843 4 : CPASSERT(.NOT. do_kpoints)
844 4 : t0_state => wfi_get_snapshot(wf_history, wf_index=2)
845 4 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
846 4 : CPASSERT(ASSOCIATED(t0_state))
847 4 : CPASSERT(ASSOCIATED(t1_state))
848 4 : CPASSERT(ASSOCIATED(t0_state%wf))
849 4 : CPASSERT(ASSOCIATED(t1_state%wf))
850 4 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
851 4 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
852 :
853 4 : my_orthogonal_wf = .TRUE.
854 4 : t0 = 0.0_dp
855 4 : t1 = t1_state%dt
856 4 : t2 = t1 + dt
857 4 : CALL qs_rho_get(rho, rho_ao=rho_ao)
858 8 : DO ispin = 1, SIZE(mos)
859 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, &
860 4 : nmo=nmo)
861 : CALL cp_fm_scale_and_add(alpha=0.0_dp, &
862 : matrix_a=mo_coeff, &
863 : matrix_b=t1_state%wf(ispin), &
864 4 : beta=(t2 - t0)/(t1 - t0))
865 : ! this copy should be unnecessary
866 : CALL cp_fm_scale_and_add(alpha=1.0_dp, &
867 : matrix_a=mo_coeff, &
868 4 : beta=(t1 - t2)/(t1 - t0), matrix_b=t0_state%wf(ispin))
869 : CALL reorthogonalize_vectors(qs_env, &
870 : v_matrix=mo_coeff, &
871 4 : n_col=nmo)
872 : CALL calculate_density_matrix(mo_set=mos(ispin), &
873 12 : density_matrix=rho_ao(ispin)%matrix)
874 : END DO
875 4 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
876 :
877 : CALL qs_ks_did_change(qs_env%ks_env, &
878 4 : rho_changed=.TRUE.)
879 : CASE (wfi_linear_p_method_nr)
880 6 : t0_state => wfi_get_snapshot(wf_history, wf_index=2)
881 6 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
882 6 : CPASSERT(ASSOCIATED(t0_state))
883 6 : CPASSERT(ASSOCIATED(t1_state))
884 6 : IF (do_kpoints) THEN
885 2 : CPASSERT(ASSOCIATED(t0_state%rho_ao_kp))
886 2 : CPASSERT(ASSOCIATED(t1_state%rho_ao_kp))
887 : ELSE
888 4 : CPASSERT(ASSOCIATED(t0_state%rho_ao))
889 4 : CPASSERT(ASSOCIATED(t1_state%rho_ao))
890 : END IF
891 6 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
892 6 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
893 :
894 6 : t0 = 0.0_dp
895 6 : t1 = t1_state%dt
896 6 : t2 = t1 + dt
897 6 : IF (do_kpoints) THEN
898 2 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
899 4 : DO ispin = 1, SIZE(rho_ao_kp, 1)
900 528 : DO img = 1, SIZE(rho_ao_kp, 2)
901 524 : IF (img > SIZE(t0_state%rho_ao_kp, 2) .OR. &
902 2 : img > SIZE(t1_state%rho_ao_kp, 2)) THEN
903 22 : CPWARN("Change in cell neighborlist: might affect quality of initial guess")
904 : ELSE
905 : CALL dbcsr_add(rho_ao_kp(ispin, img)%matrix, t1_state%rho_ao_kp(ispin, img)%matrix, &
906 502 : alpha_scalar=0.0_dp, beta_scalar=(t2 - t0)/(t1 - t0)) ! this copy should be unnecessary
907 : CALL dbcsr_add(rho_ao_kp(ispin, img)%matrix, t0_state%rho_ao_kp(ispin, img)%matrix, &
908 502 : alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0))
909 : END IF
910 : END DO
911 : END DO
912 : ELSE
913 4 : CALL qs_rho_get(rho, rho_ao=rho_ao)
914 8 : DO ispin = 1, SIZE(rho_ao)
915 : CALL dbcsr_add(rho_ao(ispin)%matrix, t1_state%rho_ao(ispin)%matrix, &
916 4 : alpha_scalar=0.0_dp, beta_scalar=(t2 - t0)/(t1 - t0)) ! this copy should be unnecessary
917 : CALL dbcsr_add(rho_ao(ispin)%matrix, t0_state%rho_ao(ispin)%matrix, &
918 8 : alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0))
919 : END DO
920 : END IF
921 6 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
922 6 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
923 : CASE (wfi_linear_ps_method_nr)
924 : ! wf not calculated, extract with PSC renormalized?
925 : ! use wf_linear?
926 12 : CPASSERT(.NOT. do_kpoints)
927 12 : t0_state => wfi_get_snapshot(wf_history, wf_index=2)
928 12 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
929 12 : CPASSERT(ASSOCIATED(t0_state))
930 12 : CPASSERT(ASSOCIATED(t1_state))
931 12 : CPASSERT(ASSOCIATED(t0_state%wf))
932 12 : CPASSERT(ASSOCIATED(t1_state%wf))
933 12 : IF (wf_history%store_overlap) THEN
934 4 : CPASSERT(ASSOCIATED(t0_state%overlap))
935 4 : CPASSERT(ASSOCIATED(t1_state%overlap))
936 : END IF
937 12 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
938 12 : IF (nvec >= wf_history%memory_depth) THEN
939 12 : IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. (qs_env%scf_control%eps_scf_hist /= 0)) THEN
940 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
941 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
942 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
943 12 : ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
944 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
945 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
946 12 : ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
947 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
948 : END IF
949 : END IF
950 :
951 12 : my_orthogonal_wf = .TRUE.
952 : ! use PS_2=2 PS_1-PS_0
953 : ! C_2 comes from using PS_2 as a projector acting on C_1
954 12 : CALL qs_rho_get(rho, rho_ao=rho_ao)
955 24 : DO ispin = 1, SIZE(mos)
956 12 : NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
957 12 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
958 : CALL cp_fm_get_info(mo_coeff, nrow_global=n, ncol_global=k, &
959 12 : matrix_struct=matrix_struct)
960 : CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, &
961 12 : nrow_global=k, ncol_global=k)
962 12 : CALL cp_fm_create(csc, matrix_struct_new)
963 12 : CALL cp_fm_struct_release(matrix_struct_new)
964 :
965 12 : IF (use_overlap) THEN
966 4 : CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), mo_coeff, k)
967 4 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), mo_coeff, 0.0_dp, csc)
968 : ELSE
969 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
970 8 : t1_state%wf(ispin), 0.0_dp, csc)
971 : END IF
972 12 : CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, mo_coeff)
973 12 : CALL cp_fm_release(csc)
974 12 : CALL cp_fm_scale_and_add(-1.0_dp, mo_coeff, 2.0_dp, t1_state%wf(ispin))
975 : CALL reorthogonalize_vectors(qs_env, &
976 : v_matrix=mo_coeff, &
977 12 : n_col=k)
978 : CALL calculate_density_matrix(mo_set=mos(ispin), &
979 48 : density_matrix=rho_ao(ispin)%matrix)
980 : END DO
981 12 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
982 12 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
983 :
984 : CASE (wfi_ps_method_nr)
985 : ! figure out the actual number of vectors to use in the extrapolation:
986 376 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
987 376 : CPASSERT(nvec > 0)
988 376 : IF (nvec >= wf_history%memory_depth) THEN
989 178 : IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. (qs_env%scf_control%eps_scf_hist /= 0)) THEN
990 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
991 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
992 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
993 178 : ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
994 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
995 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
996 178 : ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
997 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
998 : END IF
999 : END IF
1000 :
1001 376 : IF (do_kpoints) THEN
1002 4 : CALL wfi_extrapolate_ps_aspc_kp(wf_history, qs_env, nvec, io_unit, print_level)
1003 4 : my_orthogonal_wf = .TRUE.
1004 : ELSE
1005 372 : my_orthogonal_wf = .TRUE.
1006 822 : DO ispin = 1, SIZE(mos)
1007 450 : NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
1008 450 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1009 : CALL cp_fm_get_info(mo_coeff, nrow_global=n, ncol_global=k, &
1010 450 : matrix_struct=matrix_struct)
1011 450 : CALL cp_fm_create(fm_tmp, matrix_struct)
1012 : CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, &
1013 450 : nrow_global=k, ncol_global=k)
1014 450 : CALL cp_fm_create(csc, matrix_struct_new)
1015 450 : CALL cp_fm_struct_release(matrix_struct_new)
1016 : ! first the most recent
1017 450 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
1018 450 : CALL cp_fm_to_fm(t1_state%wf(ispin), mo_coeff)
1019 450 : alpha = nvec
1020 450 : CALL cp_fm_scale(alpha, mo_coeff)
1021 450 : CALL qs_rho_get(rho, rho_ao=rho_ao)
1022 962 : DO i = 2, nvec
1023 512 : t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1024 512 : IF (use_overlap) THEN
1025 474 : CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
1026 474 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
1027 : ELSE
1028 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
1029 38 : t1_state%wf(ispin), 0.0_dp, csc)
1030 : END IF
1031 512 : CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
1032 512 : alpha = -1.0_dp*alpha*REAL(nvec - i + 1, dp)/REAL(i, dp)
1033 962 : CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp)
1034 : END DO
1035 :
1036 450 : CALL cp_fm_release(csc)
1037 450 : CALL cp_fm_release(fm_tmp)
1038 : CALL reorthogonalize_vectors(qs_env, &
1039 : v_matrix=mo_coeff, &
1040 450 : n_col=k)
1041 : CALL calculate_density_matrix(mo_set=mos(ispin), &
1042 1722 : density_matrix=rho_ao(ispin)%matrix)
1043 : END DO
1044 372 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1045 372 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
1046 : END IF
1047 :
1048 : CASE (wfi_aspc_nr)
1049 13058 : CALL cite_reference(Kolafa2004)
1050 13058 : CALL cite_reference(Kuhne2007)
1051 : ! figure out the actual number of vectors to use in the extrapolation:
1052 13058 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
1053 13058 : CPASSERT(nvec > 0)
1054 13058 : IF (nvec >= wf_history%memory_depth) THEN
1055 8356 : IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. &
1056 : (qs_env%scf_control%eps_scf_hist /= 0)) THEN
1057 18 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1058 18 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1059 18 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1060 8338 : ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
1061 62 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1062 62 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1063 8276 : ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
1064 8 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1065 : END IF
1066 : END IF
1067 :
1068 13058 : IF (do_kpoints) THEN
1069 358 : CALL wfi_extrapolate_ps_aspc_kp(wf_history, qs_env, nvec, io_unit, print_level)
1070 358 : my_orthogonal_wf = .TRUE.
1071 : ELSE
1072 12700 : my_orthogonal_wf = .TRUE.
1073 12700 : CALL qs_rho_get(rho, rho_ao=rho_ao)
1074 25971 : DO ispin = 1, SIZE(mos)
1075 13271 : NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
1076 13271 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1077 : CALL cp_fm_get_info(mo_coeff, &
1078 : nrow_global=n, &
1079 : ncol_global=k, &
1080 13271 : matrix_struct=matrix_struct)
1081 13271 : CALL cp_fm_create(fm_tmp, matrix_struct, set_zero=.TRUE.)
1082 : CALL cp_fm_struct_create(matrix_struct_new, &
1083 : template_fmstruct=matrix_struct, &
1084 : nrow_global=k, &
1085 13271 : ncol_global=k)
1086 13271 : CALL cp_fm_create(csc, matrix_struct_new, set_zero=.TRUE.)
1087 13271 : CALL cp_fm_struct_release(matrix_struct_new)
1088 : ! first the most recent
1089 : t1_state => wfi_get_snapshot(wf_history, &
1090 13271 : wf_index=1)
1091 13271 : CALL cp_fm_to_fm(t1_state%wf(ispin), mo_coeff)
1092 13271 : alpha = REAL(4*nvec - 2, KIND=dp)/REAL(nvec + 1, KIND=dp)
1093 13271 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1094 : WRITE (UNIT=io_unit, FMT="(/,T2,A,/,/,T3,A,I0,/,/,T3,A2,I0,A4,F10.6)") &
1095 2500 : "Parameters for the always stable predictor-corrector (ASPC) method:", &
1096 2500 : "ASPC order: ", MAX(nvec - 2, 0), &
1097 5000 : "B(", 1, ") = ", alpha
1098 : END IF
1099 13271 : CALL cp_fm_scale(alpha, mo_coeff)
1100 :
1101 51063 : DO i = 2, nvec
1102 37792 : t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1103 37792 : IF (use_overlap) THEN
1104 26252 : CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
1105 26252 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
1106 : ELSE
1107 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
1108 11540 : t1_state%wf(ispin), 0.0_dp, csc)
1109 : END IF
1110 37792 : CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
1111 : alpha = (-1.0_dp)**(i + 1)*REAL(i, KIND=dp)* &
1112 37792 : binomial(2*nvec, nvec - i)/binomial(2*nvec - 2, nvec - 1)
1113 37792 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1114 : WRITE (UNIT=io_unit, FMT="(T3,A2,I0,A4,F10.6)") &
1115 7202 : "B(", i, ") = ", alpha
1116 : END IF
1117 51063 : CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp)
1118 : END DO
1119 13271 : CALL cp_fm_release(csc)
1120 13271 : CALL cp_fm_release(fm_tmp)
1121 : CALL reorthogonalize_vectors(qs_env, &
1122 : v_matrix=mo_coeff, &
1123 13271 : n_col=k)
1124 : CALL calculate_density_matrix(mo_set=mos(ispin), &
1125 39242 : density_matrix=rho_ao(ispin)%matrix)
1126 : END DO
1127 12700 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1128 12700 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
1129 : END IF ! do_kpoints
1130 :
1131 : CASE default
1132 : CALL cp_abort(__LOCATION__, &
1133 : "Unknown interpolation method: "// &
1134 21169 : TRIM(ADJUSTL(cp_to_string(wf_history%interpolation_method_nr))))
1135 : END SELECT
1136 21169 : IF (PRESENT(orthogonal_wf)) orthogonal_wf = my_orthogonal_wf
1137 : CALL cp_print_key_finished_output(io_unit, logger, qs_env%input, &
1138 21169 : "DFT%SCF%PRINT%PROGRAM_RUN_INFO")
1139 21169 : CALL timestop(handle)
1140 21169 : END SUBROUTINE wfi_extrapolate
1141 :
1142 : ! **************************************************************************************************
1143 : !> \brief Reorthogonalizes the wavefunctions from the previous step for k-points
1144 : !> using the current S(k) metric and rebuilds the density matrix.
1145 : !> \param qs_env The QS environment
1146 : !> \param io_unit output unit
1147 : !> \param print_level print level
1148 : ! **************************************************************************************************
1149 368 : SUBROUTINE wfi_use_prev_wf_kp(qs_env, io_unit, print_level)
1150 : TYPE(qs_environment_type), POINTER :: qs_env
1151 : INTEGER, INTENT(IN) :: io_unit, print_level
1152 :
1153 : CHARACTER(len=*), PARAMETER :: routineN = 'wfi_use_prev_wf_kp'
1154 :
1155 368 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: col_scaling
1156 : INTEGER :: chol_info, handle, igroup, ik, ikp, &
1157 : indx, ispin, j, kplocal, nao, nkp, &
1158 : nkp_groups, nmo, nspin
1159 : INTEGER, DIMENSION(2) :: kp_range
1160 368 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
1161 368 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1162 : LOGICAL :: my_kpgrp, use_real_wfn
1163 : REAL(KIND=dp) :: eval_thresh
1164 368 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1165 368 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1166 368 : TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
1167 : TYPE(cp_cfm_type) :: cfm_evecs, cfm_mhalf, cfm_nao_nmo_work, &
1168 : cmos_new, csc_cfm
1169 368 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: csmat_cur
1170 368 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools_kp
1171 : TYPE(cp_fm_struct_type), POINTER :: ao_ao_struct, nmo_nmo_struct
1172 : TYPE(cp_fm_type) :: fmdummy, fmlocal
1173 : TYPE(cp_fm_type), POINTER :: imos, mo_coeff, rmos
1174 368 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, rho_ao_kp
1175 : TYPE(dbcsr_type), POINTER :: cmatrix_db, rmatrix, tmpmat
1176 : TYPE(dft_control_type), POINTER :: dft_control
1177 : TYPE(kpoint_env_type), POINTER :: kp
1178 : TYPE(kpoint_type), POINTER :: kpoints
1179 : TYPE(mp_para_env_type), POINTER :: para_env
1180 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1181 368 : POINTER :: sab_nl
1182 : TYPE(qs_matrix_pools_type), POINTER :: mpools_kp
1183 : TYPE(qs_rho_type), POINTER :: rho
1184 : TYPE(qs_scf_env_type), POINTER :: scf_env
1185 : TYPE(scf_control_type), POINTER :: scf_control
1186 :
1187 368 : CALL timeset(routineN, handle)
1188 :
1189 368 : NULLIFY (kpoints, dft_control, matrix_s_kp, scf_env, scf_control, rho, sab_nl, kp, mo_coeff, rmos, imos)
1190 :
1191 : CALL get_qs_env(qs_env, kpoints=kpoints, dft_control=dft_control, &
1192 : matrix_s_kp=matrix_s_kp, scf_env=scf_env, &
1193 368 : scf_control=scf_control, rho=rho)
1194 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn, &
1195 : kp_range=kp_range, nkp_groups=nkp_groups, kp_dist=kp_dist, &
1196 368 : sab_nl=sab_nl, cell_to_index=cell_to_index)
1197 368 : kplocal = kp_range(2) - kp_range(1) + 1
1198 :
1199 368 : IF (use_real_wfn) THEN
1200 0 : CALL timestop(handle)
1201 0 : RETURN
1202 : END IF
1203 :
1204 368 : kp => kpoints%kp_env(1)%kpoint_env
1205 368 : nspin = SIZE(kp%mos, 2)
1206 368 : CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo, mo_coeff=mo_coeff)
1207 :
1208 368 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1209 : WRITE (UNIT=io_unit, FMT="(/,T2,A)") &
1210 0 : "Using previous wavefunctions as initial guess for k-points (with reorthogonalization)"
1211 : END IF
1212 :
1213 : ! Allocate dbcsr work matrices
1214 368 : ALLOCATE (rmatrix, cmatrix_db, tmpmat)
1215 368 : CALL dbcsr_create(rmatrix, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
1216 368 : CALL dbcsr_create(cmatrix_db, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
1217 368 : CALL dbcsr_create(tmpmat, template=matrix_s_kp(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1218 368 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
1219 368 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix_db, sab_nl)
1220 :
1221 368 : CALL get_kpoint_info(kpoints, mpools=mpools_kp)
1222 368 : CALL mpools_get(mpools_kp, ao_ao_fm_pools=ao_ao_fm_pools_kp)
1223 368 : CALL fm_pool_create_fm(ao_ao_fm_pools_kp(1)%pool, fmlocal)
1224 368 : CALL cp_fm_get_info(fmlocal, matrix_struct=ao_ao_struct)
1225 :
1226 368 : CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct)
1227 368 : CALL cp_cfm_create(cfm_nao_nmo_work, mo_coeff%matrix_struct)
1228 :
1229 368 : NULLIFY (nmo_nmo_struct)
1230 : CALL cp_fm_struct_create(nmo_nmo_struct, template_fmstruct=mo_coeff%matrix_struct, &
1231 368 : nrow_global=nmo, ncol_global=nmo)
1232 368 : CALL cp_cfm_create(csc_cfm, nmo_nmo_struct)
1233 368 : CALL cp_fm_struct_release(nmo_nmo_struct)
1234 :
1235 368 : para_env => kpoints%blacs_env_all%para_env
1236 7456 : ALLOCATE (info(kplocal*nkp_groups, 2))
1237 :
1238 2104 : ALLOCATE (csmat_cur(kplocal))
1239 1368 : DO ikp = 1, kplocal
1240 1368 : CALL cp_cfm_create(csmat_cur(ikp), ao_ao_struct)
1241 : END DO
1242 :
1243 : ! Phase A: Fourier Transform S(R) -> S(k)
1244 : indx = 0
1245 1368 : DO ikp = 1, kplocal
1246 3072 : DO igroup = 1, nkp_groups
1247 1704 : ik = kp_dist(1, igroup) + ikp - 1
1248 1704 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
1249 1704 : indx = indx + 1
1250 :
1251 1704 : CALL dbcsr_set(rmatrix, 0.0_dp)
1252 1704 : CALL dbcsr_set(cmatrix_db, 0.0_dp)
1253 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix_db, rsmat=matrix_s_kp, &
1254 1704 : ispin=1, xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
1255 1704 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
1256 1704 : CALL copy_dbcsr_to_fm(tmpmat, scf_env%scf_work1(1))
1257 1704 : CALL dbcsr_desymmetrize(cmatrix_db, tmpmat)
1258 1704 : CALL copy_dbcsr_to_fm(tmpmat, scf_env%scf_work1(2))
1259 :
1260 2704 : IF (my_kpgrp) THEN
1261 1000 : CALL cp_fm_start_copy_general(scf_env%scf_work1(1), fmlocal, para_env, info(indx, 1))
1262 1000 : CALL cp_fm_start_copy_general(scf_env%scf_work1(2), fmlocal, para_env, info(indx, 2))
1263 : ELSE
1264 704 : CALL cp_fm_start_copy_general(scf_env%scf_work1(1), fmdummy, para_env, info(indx, 1))
1265 704 : CALL cp_fm_start_copy_general(scf_env%scf_work1(2), fmdummy, para_env, info(indx, 2))
1266 : END IF
1267 : END DO
1268 : END DO
1269 :
1270 : ! Finish Communication
1271 : indx = 0
1272 1368 : DO ikp = 1, kplocal
1273 3072 : DO igroup = 1, nkp_groups
1274 1704 : ik = kp_dist(1, igroup) + ikp - 1
1275 1704 : my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2))
1276 704 : indx = indx + 1
1277 1000 : IF (my_kpgrp) THEN
1278 1000 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
1279 1000 : CALL cp_cfm_scale_and_add_fm(z_zero, csmat_cur(ikp), z_one, fmlocal)
1280 1000 : CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
1281 1000 : CALL cp_cfm_scale_and_add_fm(z_one, csmat_cur(ikp), gaussi, fmlocal)
1282 : END IF
1283 : END DO
1284 : END DO
1285 :
1286 2072 : DO indx = 1, kplocal*nkp_groups
1287 1704 : CALL cp_fm_cleanup_copy_general(info(indx, 1))
1288 2072 : CALL cp_fm_cleanup_copy_general(info(indx, 2))
1289 : END DO
1290 :
1291 : ! Phase B: Orthogonalize current MOs w.r.t. new S(k)
1292 1104 : ALLOCATE (eigenvalues(nmo))
1293 368 : eval_thresh = 1.0E-12_dp
1294 :
1295 1368 : DO ikp = 1, kplocal
1296 1000 : kp => kpoints%kp_env(ikp)%kpoint_env
1297 2560 : DO ispin = 1, nspin
1298 1192 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
1299 1192 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1300 1192 : CALL cp_fm_to_cfm(rmos, imos, cmos_new)
1301 :
1302 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nao, z_one, &
1303 1192 : csmat_cur(ikp), cmos_new, z_zero, cfm_nao_nmo_work)
1304 : CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, &
1305 1192 : cmos_new, cfm_nao_nmo_work, z_zero, csc_cfm)
1306 :
1307 1192 : CALL cp_cfm_cholesky_decompose(csc_cfm, info_out=chol_info)
1308 1192 : IF (chol_info == 0) THEN
1309 1192 : CALL cp_cfm_triangular_multiply(csc_cfm, cmos_new, side='R', invert_tr=.TRUE., uplo_tr='U')
1310 : ELSE
1311 0 : CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, cmos_new, cfm_nao_nmo_work, z_zero, csc_cfm)
1312 0 : CALL cp_cfm_create(cfm_evecs, csc_cfm%matrix_struct)
1313 0 : CALL cp_cfm_create(cfm_mhalf, csc_cfm%matrix_struct)
1314 0 : CALL cp_cfm_heevd(csc_cfm, cfm_evecs, eigenvalues)
1315 0 : CALL cp_cfm_to_cfm(cfm_evecs, cfm_mhalf)
1316 0 : ALLOCATE (col_scaling(nmo))
1317 0 : DO j = 1, nmo
1318 0 : IF (eigenvalues(j) > eval_thresh) THEN
1319 0 : col_scaling(j) = CMPLX(1.0_dp/SQRT(eigenvalues(j)), 0.0_dp, KIND=dp)
1320 : ELSE
1321 0 : col_scaling(j) = z_zero
1322 : END IF
1323 : END DO
1324 0 : CALL cp_cfm_column_scale(cfm_mhalf, col_scaling)
1325 0 : DEALLOCATE (col_scaling)
1326 0 : CALL cp_cfm_gemm('N', 'C', nmo, nmo, nmo, z_one, cfm_mhalf, cfm_evecs, z_zero, csc_cfm)
1327 0 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nmo, z_one, cmos_new, csc_cfm, z_zero, cfm_nao_nmo_work)
1328 0 : CALL cp_cfm_to_cfm(cfm_nao_nmo_work, cmos_new)
1329 0 : CALL cp_cfm_release(cfm_evecs)
1330 0 : CALL cp_cfm_release(cfm_mhalf)
1331 : END IF
1332 3384 : CALL cp_cfm_to_fm(cmos_new, rmos, imos)
1333 : END DO
1334 : END DO
1335 368 : DEALLOCATE (eigenvalues)
1336 :
1337 : ! Phase C: Rebuild Density Matrix P(R)
1338 368 : CALL kpoint_set_mo_occupation(kpoints, scf_control%smear)
1339 368 : CALL kpoint_density_matrices(kpoints)
1340 368 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1341 : CALL kpoint_density_transform(kpoints, rho_ao_kp, .FALSE., &
1342 368 : matrix_s_kp(1, 1)%matrix, sab_nl, scf_env%scf_work1)
1343 368 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
1344 368 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
1345 :
1346 : ! Cleanup
1347 1368 : DO ikp = 1, kplocal
1348 1368 : CALL cp_cfm_release(csmat_cur(ikp))
1349 : END DO
1350 368 : DEALLOCATE (csmat_cur)
1351 3776 : DEALLOCATE (info)
1352 368 : CALL cp_cfm_release(cmos_new)
1353 368 : CALL cp_cfm_release(cfm_nao_nmo_work)
1354 368 : CALL cp_cfm_release(csc_cfm)
1355 368 : CALL fm_pool_give_back_fm(ao_ao_fm_pools_kp(1)%pool, fmlocal)
1356 368 : CALL dbcsr_deallocate_matrix(rmatrix)
1357 368 : CALL dbcsr_deallocate_matrix(cmatrix_db)
1358 368 : CALL dbcsr_deallocate_matrix(tmpmat)
1359 :
1360 368 : CALL timestop(handle)
1361 2208 : END SUBROUTINE wfi_use_prev_wf_kp
1362 :
1363 : ! **************************************************************************************************
1364 : !> \brief Performs PS/ASPC wavefunction extrapolation for k-point calculations.
1365 : !> Applies PS/ASPC coefficients to complex MO coefficients at each k-point,
1366 : !> with subspace alignment via historical overlap matrices.
1367 : !> Delegates final orthogonalization and density building to wfi_use_prev_wf_kp.
1368 : !> \param wf_history wavefunction history buffer
1369 : !> \param qs_env QS environment
1370 : !> \param nvec number of history snapshots to use
1371 : !> \param io_unit output unit for logging
1372 : !> \param print_level current print level
1373 : ! **************************************************************************************************
1374 724 : SUBROUTINE wfi_extrapolate_ps_aspc_kp(wf_history, qs_env, nvec, io_unit, print_level)
1375 : TYPE(qs_wf_history_type), POINTER :: wf_history
1376 : TYPE(qs_environment_type), POINTER :: qs_env
1377 : INTEGER, INTENT(IN) :: nvec, io_unit, print_level
1378 :
1379 : CHARACTER(len=*), PARAMETER :: routineN = 'wfi_extrapolate_ps_aspc_kp'
1380 :
1381 : INTEGER :: handle, i, ikp, ispin, kplocal, &
1382 : method_nr, nao, nmo, nspin
1383 : INTEGER, DIMENSION(2) :: kp_range
1384 : LOGICAL :: use_real_wfn
1385 : REAL(KIND=dp) :: alpha_coeff
1386 : TYPE(cp_cfm_type) :: cfm_nao_nmo_work, cmos_1, cmos_i, &
1387 : cmos_new, csc_cfm
1388 : TYPE(cp_fm_struct_type), POINTER :: nmo_nmo_struct
1389 : TYPE(cp_fm_type), POINTER :: imos, mo_coeff, rmos
1390 : TYPE(kpoint_env_type), POINTER :: kp
1391 : TYPE(kpoint_type), POINTER :: kpoints
1392 : TYPE(qs_wf_snapshot_type), POINTER :: t0_state, t1_state
1393 :
1394 362 : method_nr = wf_history%interpolation_method_nr
1395 :
1396 362 : CALL timeset(routineN, handle)
1397 362 : NULLIFY (kpoints, kp, mo_coeff, rmos, imos, t0_state, t1_state, nmo_nmo_struct)
1398 :
1399 362 : CALL get_qs_env(qs_env, kpoints=kpoints)
1400 362 : CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn, kp_range=kp_range)
1401 362 : kplocal = kp_range(2) - kp_range(1) + 1
1402 :
1403 362 : IF (use_real_wfn) THEN
1404 0 : IF (method_nr == wfi_aspc_nr) THEN
1405 : CALL cp_warn(__LOCATION__, "ASPC with k-points requires complex wavefunctions; "// &
1406 0 : "falling back to USE_PREV_WF.")
1407 : ELSE
1408 : CALL cp_warn(__LOCATION__, "PS with k-points requires complex wavefunctions; "// &
1409 0 : "falling back to USE_PREV_WF.")
1410 : END IF
1411 0 : CALL wfi_use_prev_wf_kp(qs_env, io_unit, print_level)
1412 0 : CALL timestop(handle)
1413 0 : RETURN
1414 : END IF
1415 :
1416 362 : kp => kpoints%kp_env(1)%kpoint_env
1417 362 : nspin = SIZE(kp%mos, 2)
1418 362 : CALL get_mo_set(kp%mos(1, 1), nao=nao, nmo=nmo, mo_coeff=mo_coeff)
1419 :
1420 362 : IF (method_nr == wfi_aspc_nr) THEN
1421 358 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1422 : WRITE (UNIT=io_unit, FMT="(/,T2,A,/,T3,A,I0)") &
1423 7 : "Parameters for the always stable predictor-corrector (ASPC) method:", &
1424 14 : "ASPC order: ", MAX(nvec - 2, 0)
1425 : END IF
1426 : END IF
1427 :
1428 16 : IF (method_nr == wfi_aspc_nr) THEN
1429 358 : CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct, set_zero=.TRUE.)
1430 358 : CALL cp_cfm_create(cmos_1, mo_coeff%matrix_struct, set_zero=.TRUE.)
1431 358 : CALL cp_cfm_create(cmos_i, mo_coeff%matrix_struct, set_zero=.TRUE.)
1432 358 : CALL cp_cfm_create(cfm_nao_nmo_work, mo_coeff%matrix_struct, set_zero=.TRUE.)
1433 : ELSE
1434 4 : CALL cp_cfm_create(cmos_new, mo_coeff%matrix_struct)
1435 4 : CALL cp_cfm_create(cmos_1, mo_coeff%matrix_struct)
1436 4 : CALL cp_cfm_create(cmos_i, mo_coeff%matrix_struct)
1437 4 : CALL cp_cfm_create(cfm_nao_nmo_work, mo_coeff%matrix_struct)
1438 : END IF
1439 :
1440 : CALL cp_fm_struct_create(nmo_nmo_struct, template_fmstruct=mo_coeff%matrix_struct, &
1441 362 : nrow_global=nmo, ncol_global=nmo)
1442 362 : IF (method_nr == wfi_aspc_nr) THEN
1443 358 : CALL cp_cfm_create(csc_cfm, nmo_nmo_struct, set_zero=.TRUE.)
1444 : ELSE
1445 4 : CALL cp_cfm_create(csc_cfm, nmo_nmo_struct)
1446 : END IF
1447 362 : CALL cp_fm_struct_release(nmo_nmo_struct)
1448 :
1449 : ! Phase 1: Initialize C_new(k) = B(1) * C_1(k)
1450 362 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
1451 362 : IF (method_nr == wfi_aspc_nr) THEN
1452 358 : alpha_coeff = REAL(4*nvec - 2, KIND=dp)/REAL(nvec + 1, KIND=dp)
1453 358 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1454 7 : WRITE (UNIT=io_unit, FMT="(T3,A2,I0,A4,F10.6)") "B(", 1, ") = ", alpha_coeff
1455 : END IF
1456 : ELSE
1457 4 : alpha_coeff = nvec
1458 : END IF
1459 :
1460 1338 : DO ikp = 1, kplocal
1461 976 : kp => kpoints%kp_env(ikp)%kpoint_env
1462 2506 : DO ispin = 1, nspin
1463 1168 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
1464 1168 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1465 1168 : CALL cp_fm_to_fm(t1_state%wf_kp(ikp, 1, ispin), rmos)
1466 1168 : CALL cp_fm_to_fm(t1_state%wf_kp(ikp, 2, ispin), imos)
1467 1168 : CALL cp_fm_scale(alpha_coeff, rmos)
1468 2144 : CALL cp_fm_scale(alpha_coeff, imos)
1469 : END DO
1470 : END DO
1471 :
1472 : ! Phase 2: Accumulate historical snapshots C_new += B(i) * C_proj(k)
1473 1544 : DO i = 2, nvec
1474 1182 : t0_state => wfi_get_snapshot(wf_history, wf_index=i)
1475 1182 : IF (method_nr == wfi_aspc_nr) THEN
1476 : alpha_coeff = (-1.0_dp)**(i + 1)*REAL(i, KIND=dp)* &
1477 1180 : binomial(2*nvec, nvec - i)/binomial(2*nvec - 2, nvec - 1)
1478 1180 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
1479 5 : WRITE (UNIT=io_unit, FMT="(T3,A2,I0,A4,F10.6)") "B(", i, ") = ", alpha_coeff
1480 : END IF
1481 : ELSE
1482 2 : alpha_coeff = -1.0_dp*alpha_coeff*REAL(nvec - i + 1, dp)/REAL(i, dp)
1483 : END IF
1484 :
1485 3816 : DO ikp = 1, kplocal
1486 2272 : kp => kpoints%kp_env(ikp)%kpoint_env
1487 5822 : DO ispin = 1, nspin
1488 2368 : CALL cp_fm_to_cfm(t1_state%wf_kp(ikp, 1, ispin), t1_state%wf_kp(ikp, 2, ispin), cmos_1)
1489 2368 : CALL cp_fm_to_cfm(t0_state%wf_kp(ikp, 1, ispin), t0_state%wf_kp(ikp, 2, ispin), cmos_i)
1490 :
1491 : ! Subspace projection: C_proj = C_i * (C_i^dag S_i C_1)
1492 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nao, z_one, &
1493 2368 : t0_state%overlap_cfm_kp(ikp), cmos_1, z_zero, cfm_nao_nmo_work)
1494 : CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, &
1495 2368 : cmos_i, cfm_nao_nmo_work, z_zero, csc_cfm)
1496 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nmo, z_one, &
1497 2368 : cmos_i, csc_cfm, z_zero, cfm_nao_nmo_work)
1498 :
1499 2368 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos)
1500 2368 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1501 2368 : CALL cp_fm_to_cfm(rmos, imos, cmos_new)
1502 2368 : CALL cp_cfm_scale_and_add(z_one, cmos_new, CMPLX(alpha_coeff, 0.0_dp, KIND=dp), cfm_nao_nmo_work)
1503 4640 : CALL cp_cfm_to_fm(cmos_new, rmos, imos)
1504 : END DO
1505 : END DO
1506 : END DO
1507 :
1508 362 : CALL cp_cfm_release(cmos_new)
1509 362 : CALL cp_cfm_release(cmos_1)
1510 362 : CALL cp_cfm_release(cmos_i)
1511 362 : CALL cp_cfm_release(cfm_nao_nmo_work)
1512 362 : CALL cp_cfm_release(csc_cfm)
1513 :
1514 : ! Phase 3: Delegate Orthogonalization and Density Building
1515 : ! We pass io_unit = 0 to suppress the redundant "Using previous..." print
1516 : ! inside wfi_use_prev_wf_kp, keeping the ASPC log output perfectly clean.
1517 362 : CALL wfi_use_prev_wf_kp(qs_env, 0, print_level)
1518 :
1519 362 : CALL timestop(handle)
1520 :
1521 362 : END SUBROUTINE wfi_extrapolate_ps_aspc_kp
1522 :
1523 : ! **************************************************************************************************
1524 : !> \brief Decides if scf control variables has to changed due
1525 : !> to using a WF extrapolation.
1526 : !> \param qs_env The QS environment
1527 : !> \param nvec ...
1528 : !> \par History
1529 : !> 11.2006 created [TdK]
1530 : !> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch)
1531 : ! **************************************************************************************************
1532 7723 : ELEMENTAL SUBROUTINE wfi_set_history_variables(qs_env, nvec)
1533 : TYPE(qs_environment_type), INTENT(INOUT) :: qs_env
1534 : INTEGER, INTENT(IN) :: nvec
1535 :
1536 7723 : IF (nvec >= qs_env%wf_history%memory_depth) THEN
1537 1289 : IF ((qs_env%scf_control%max_scf_hist /= 0) .AND. (qs_env%scf_control%eps_scf_hist /= 0)) THEN
1538 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1539 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1540 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1541 1289 : ELSE IF (qs_env%scf_control%max_scf_hist /= 0) THEN
1542 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
1543 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
1544 1289 : ELSE IF (qs_env%scf_control%eps_scf_hist /= 0) THEN
1545 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
1546 0 : qs_env%scf_control%outer_scf%eps_scf = qs_env%scf_control%eps_scf_hist
1547 : END IF
1548 : END IF
1549 :
1550 7723 : END SUBROUTINE wfi_set_history_variables
1551 :
1552 : ! **************************************************************************************************
1553 : !> \brief updates the snapshot buffer, taking a new snapshot
1554 : !> \param wf_history the history buffer to update
1555 : !> \param qs_env the qs_env we get the info from
1556 : !> \param dt ...
1557 : !> \par History
1558 : !> 02.2003 created [fawzi]
1559 : !> \author fawzi
1560 : ! **************************************************************************************************
1561 21173 : SUBROUTINE wfi_update(wf_history, qs_env, dt)
1562 : TYPE(qs_wf_history_type), POINTER :: wf_history
1563 : TYPE(qs_environment_type), POINTER :: qs_env
1564 : REAL(KIND=dp), INTENT(in) :: dt
1565 :
1566 21173 : CPASSERT(ASSOCIATED(wf_history))
1567 21173 : CPASSERT(wf_history%ref_count > 0)
1568 21173 : CPASSERT(ASSOCIATED(qs_env))
1569 :
1570 21173 : wf_history%snapshot_count = wf_history%snapshot_count + 1
1571 21173 : IF (wf_history%memory_depth > 0) THEN
1572 : wf_history%last_state_index = MODULO(wf_history%snapshot_count, &
1573 20386 : wf_history%memory_depth) + 1
1574 : CALL wfs_update(snapshot=wf_history%past_states &
1575 : (wf_history%last_state_index)%snapshot, wf_history=wf_history, &
1576 20386 : qs_env=qs_env, dt=dt)
1577 : END IF
1578 21173 : END SUBROUTINE wfi_update
1579 :
1580 : ! **************************************************************************************************
1581 : !> \brief reorthogonalizes the mos
1582 : !> \param qs_env the qs_env in which to orthogonalize
1583 : !> \param v_matrix the vectors to orthogonalize
1584 : !> \param n_col number of column of v to orthogonalize
1585 : !> \par History
1586 : !> 04.2003 created [fawzi]
1587 : !> \author Fawzi Mohamed
1588 : ! **************************************************************************************************
1589 28410 : SUBROUTINE reorthogonalize_vectors(qs_env, v_matrix, n_col)
1590 : TYPE(qs_environment_type), POINTER :: qs_env
1591 : TYPE(cp_fm_type), INTENT(IN) :: v_matrix
1592 : INTEGER, INTENT(in), OPTIONAL :: n_col
1593 :
1594 : CHARACTER(len=*), PARAMETER :: routineN = 'reorthogonalize_vectors'
1595 :
1596 : INTEGER :: handle, my_n_col
1597 : LOGICAL :: has_unit_metric, &
1598 : ortho_contains_cholesky, &
1599 : smearing_is_used
1600 : TYPE(cp_fm_pool_type), POINTER :: maxao_maxmo_fm_pool
1601 14205 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1602 : TYPE(dft_control_type), POINTER :: dft_control
1603 : TYPE(qs_matrix_pools_type), POINTER :: mpools
1604 : TYPE(qs_scf_env_type), POINTER :: scf_env
1605 : TYPE(scf_control_type), POINTER :: scf_control
1606 :
1607 14205 : NULLIFY (scf_env, scf_control, maxao_maxmo_fm_pool, matrix_s, mpools, dft_control)
1608 14205 : CALL timeset(routineN, handle)
1609 :
1610 14205 : CPASSERT(ASSOCIATED(qs_env))
1611 :
1612 14205 : CALL cp_fm_get_info(v_matrix, ncol_global=my_n_col)
1613 14205 : IF (PRESENT(n_col)) my_n_col = n_col
1614 : CALL get_qs_env(qs_env, mpools=mpools, &
1615 : scf_env=scf_env, &
1616 : scf_control=scf_control, &
1617 : matrix_s=matrix_s, &
1618 14205 : dft_control=dft_control)
1619 14205 : CALL mpools_get(mpools, maxao_maxmo_fm_pool=maxao_maxmo_fm_pool)
1620 14205 : IF (ASSOCIATED(scf_env)) THEN
1621 : ortho_contains_cholesky = (scf_env%method /= ot_method_nr) .AND. &
1622 : (scf_env%cholesky_method > 0) .AND. &
1623 14205 : ASSOCIATED(scf_env%ortho)
1624 : ELSE
1625 : ortho_contains_cholesky = .FALSE.
1626 : END IF
1627 :
1628 14205 : CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
1629 14205 : smearing_is_used = .FALSE.
1630 14205 : IF (dft_control%smear) THEN
1631 1148 : smearing_is_used = .TRUE.
1632 : END IF
1633 :
1634 14205 : IF (has_unit_metric) THEN
1635 3410 : CALL make_basis_simple(v_matrix, my_n_col)
1636 10795 : ELSE IF (smearing_is_used) THEN
1637 : CALL make_basis_lowdin(vmatrix=v_matrix, ncol=my_n_col, &
1638 1148 : matrix_s=matrix_s(1)%matrix)
1639 9647 : ELSE IF (ortho_contains_cholesky) THEN
1640 : CALL make_basis_cholesky(vmatrix=v_matrix, ncol=my_n_col, &
1641 6266 : ortho=scf_env%ortho)
1642 : ELSE
1643 3381 : CALL make_basis_sm(v_matrix, my_n_col, matrix_s(1)%matrix)
1644 : END IF
1645 14205 : CALL timestop(handle)
1646 14205 : END SUBROUTINE reorthogonalize_vectors
1647 :
1648 : ! **************************************************************************************************
1649 : !> \brief purges wf_history retaining only the latest snapshot
1650 : !> \param qs_env the qs env with the latest result, and that will contain
1651 : !> the purged wf_history
1652 : !> \par History
1653 : !> 05.2016 created [Nico Holmberg]
1654 : !> \author Nico Holmberg
1655 : ! **************************************************************************************************
1656 0 : SUBROUTINE wfi_purge_history(qs_env)
1657 : TYPE(qs_environment_type), POINTER :: qs_env
1658 :
1659 : CHARACTER(len=*), PARAMETER :: routineN = 'wfi_purge_history'
1660 :
1661 : INTEGER :: handle, io_unit, print_level
1662 : TYPE(cp_logger_type), POINTER :: logger
1663 : TYPE(dft_control_type), POINTER :: dft_control
1664 : TYPE(qs_wf_history_type), POINTER :: wf_history
1665 :
1666 0 : NULLIFY (dft_control, wf_history)
1667 :
1668 0 : CALL timeset(routineN, handle)
1669 0 : logger => cp_get_default_logger()
1670 0 : print_level = logger%iter_info%print_level
1671 : io_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%SCF%PRINT%PROGRAM_RUN_INFO", &
1672 0 : extension=".scfLog")
1673 :
1674 0 : CPASSERT(ASSOCIATED(qs_env))
1675 0 : CPASSERT(ASSOCIATED(qs_env%wf_history))
1676 0 : CPASSERT(qs_env%wf_history%ref_count > 0)
1677 0 : CALL get_qs_env(qs_env, dft_control=dft_control)
1678 :
1679 0 : SELECT CASE (qs_env%wf_history%interpolation_method_nr)
1680 : CASE (wfi_use_guess_method_nr, wfi_use_prev_wf_method_nr, &
1681 : wfi_use_prev_p_method_nr, wfi_use_prev_rho_r_method_nr, &
1682 : wfi_frozen_method_nr)
1683 : ! do nothing
1684 : CASE (wfi_linear_wf_method_nr, wfi_linear_p_method_nr, &
1685 : wfi_linear_ps_method_nr, wfi_ps_method_nr, &
1686 : wfi_aspc_nr)
1687 0 : IF (qs_env%wf_history%snapshot_count >= 2) THEN
1688 0 : IF (debug_this_module .AND. io_unit > 0) &
1689 0 : WRITE (io_unit, FMT="(T2,A)") "QS| Purging WFN history"
1690 : CALL wfi_create(wf_history, interpolation_method_nr= &
1691 : dft_control%qs_control%wf_interpolation_method_nr, &
1692 : extrapolation_order=dft_control%qs_control%wf_extrapolation_order, &
1693 0 : has_unit_metric=qs_env%has_unit_metric)
1694 : CALL set_qs_env(qs_env=qs_env, &
1695 0 : wf_history=wf_history)
1696 0 : CALL wfi_release(wf_history)
1697 0 : CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
1698 : END IF
1699 : CASE DEFAULT
1700 0 : CPABORT("Unknown extrapolation method.")
1701 : END SELECT
1702 0 : CALL timestop(handle)
1703 :
1704 0 : END SUBROUTINE wfi_purge_history
1705 :
1706 : END MODULE qs_wf_history_methods
|
