Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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_control_types, ONLY: dft_control_type
27 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
28 : dbcsr_copy,&
29 : dbcsr_deallocate_matrix,&
30 : dbcsr_p_type
31 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
32 : dbcsr_allocate_matrix_set,&
33 : dbcsr_deallocate_matrix_set
34 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale,&
35 : cp_fm_scale_and_add
36 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
37 : cp_fm_pool_type,&
38 : fm_pools_create_fm_vect,&
39 : fm_pools_give_back_fm_vect
40 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
41 : cp_fm_struct_release,&
42 : cp_fm_struct_type
43 : USE cp_fm_types, ONLY: cp_fm_create,&
44 : cp_fm_get_info,&
45 : cp_fm_release,&
46 : cp_fm_to_fm,&
47 : cp_fm_type
48 : USE cp_log_handling, ONLY: cp_get_default_logger,&
49 : cp_logger_type,&
50 : cp_to_string
51 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
52 : cp_print_key_unit_nr,&
53 : low_print_level
54 : USE input_constants, ONLY: &
55 : wfi_aspc_nr, wfi_frozen_method_nr, wfi_linear_p_method_nr, wfi_linear_ps_method_nr, &
56 : wfi_linear_wf_method_nr, wfi_ps_method_nr, wfi_use_guess_method_nr, &
57 : wfi_use_prev_p_method_nr, wfi_use_prev_rho_r_method_nr, wfi_use_prev_wf_method_nr
58 : USE kinds, ONLY: dp
59 : USE mathlib, ONLY: binomial
60 : USE parallel_gemm_api, ONLY: parallel_gemm
61 : USE pw_env_types, ONLY: pw_env_get,&
62 : pw_env_type
63 : USE pw_methods, ONLY: pw_copy
64 : USE pw_pool_types, ONLY: pw_pool_type
65 : USE pw_types, ONLY: pw_c1d_gs_type,&
66 : pw_r3d_rs_type
67 : USE qs_density_matrices, ONLY: calculate_density_matrix
68 : USE qs_environment_types, ONLY: get_qs_env,&
69 : qs_environment_type,&
70 : set_qs_env
71 : USE qs_ks_types, ONLY: qs_ks_did_change
72 : USE qs_matrix_pools, ONLY: mpools_get,&
73 : qs_matrix_pools_type
74 : USE qs_mo_methods, ONLY: make_basis_cholesky,&
75 : make_basis_lowdin,&
76 : make_basis_simple,&
77 : make_basis_sm
78 : USE qs_mo_types, ONLY: get_mo_set,&
79 : mo_set_type
80 : USE qs_rho_methods, ONLY: qs_rho_update_rho
81 : USE qs_rho_types, ONLY: qs_rho_get,&
82 : qs_rho_set,&
83 : qs_rho_type
84 : USE qs_scf_types, ONLY: ot_method_nr,&
85 : qs_scf_env_type
86 : USE qs_wf_history_types, ONLY: qs_wf_history_type,&
87 : qs_wf_snapshot_type,&
88 : wfi_get_snapshot,&
89 : wfi_release
90 : USE scf_control_types, ONLY: scf_control_type
91 : #include "./base/base_uses.f90"
92 :
93 : IMPLICIT NONE
94 : PRIVATE
95 :
96 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
97 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_wf_history_methods'
98 :
99 : PUBLIC :: wfi_create, wfi_update, wfi_create_for_kp, &
100 : wfi_extrapolate, wfi_get_method_label, &
101 : reorthogonalize_vectors, wfi_purge_history
102 :
103 : CONTAINS
104 :
105 : ! **************************************************************************************************
106 : !> \brief allocates and initialize a wavefunction snapshot
107 : !> \param snapshot the snapshot to create
108 : !> \par History
109 : !> 02.2003 created [fawzi]
110 : !> 02.2005 added wf_mol [MI]
111 : !> \author fawzi
112 : ! **************************************************************************************************
113 10518 : SUBROUTINE wfs_create(snapshot)
114 : TYPE(qs_wf_snapshot_type), INTENT(OUT) :: snapshot
115 :
116 : NULLIFY (snapshot%wf, snapshot%rho_r, &
117 : snapshot%rho_g, snapshot%rho_ao, snapshot%rho_ao_kp, &
118 : snapshot%overlap, snapshot%rho_frozen)
119 10518 : snapshot%dt = 1.0_dp
120 10518 : END SUBROUTINE wfs_create
121 :
122 : ! **************************************************************************************************
123 : !> \brief updates the given snapshot
124 : !> \param snapshot the snapshot to be updated
125 : !> \param wf_history the history
126 : !> \param qs_env the qs_env that should be snapshotted
127 : !> \param dt the time of the snapshot (wrt. to the previous snapshot)
128 : !> \par History
129 : !> 02.2003 created [fawzi]
130 : !> 02.2005 added kg_fm_mol_set for KG_GPW [MI]
131 : !> \author fawzi
132 : ! **************************************************************************************************
133 19078 : SUBROUTINE wfs_update(snapshot, wf_history, qs_env, dt)
134 : TYPE(qs_wf_snapshot_type), POINTER :: snapshot
135 : TYPE(qs_wf_history_type), POINTER :: wf_history
136 : TYPE(qs_environment_type), POINTER :: qs_env
137 : REAL(KIND=dp), INTENT(in), OPTIONAL :: dt
138 :
139 : CHARACTER(len=*), PARAMETER :: routineN = 'wfs_update'
140 :
141 : INTEGER :: handle, img, ispin, nimg, nspins
142 19078 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_pools
143 : TYPE(cp_fm_type), POINTER :: mo_coeff
144 19078 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
145 19078 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
146 : TYPE(dft_control_type), POINTER :: dft_control
147 19078 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
148 19078 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
149 : TYPE(pw_env_type), POINTER :: pw_env
150 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
151 19078 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
152 : TYPE(qs_rho_type), POINTER :: rho
153 :
154 19078 : CALL timeset(routineN, handle)
155 :
156 19078 : NULLIFY (pw_env, auxbas_pw_pool, ao_mo_pools, dft_control, mos, mo_coeff, &
157 19078 : rho, rho_r, rho_g, rho_ao, matrix_s)
158 : CALL get_qs_env(qs_env, pw_env=pw_env, &
159 19078 : dft_control=dft_control, rho=rho)
160 19078 : CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_pools)
161 19078 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
162 :
163 19078 : CPASSERT(ASSOCIATED(wf_history))
164 19078 : CPASSERT(ASSOCIATED(dft_control))
165 19078 : IF (.NOT. ASSOCIATED(snapshot)) THEN
166 10518 : ALLOCATE (snapshot)
167 10518 : CALL wfs_create(snapshot)
168 : END IF
169 19078 : CPASSERT(wf_history%ref_count > 0)
170 :
171 19078 : nspins = dft_control%nspins
172 19078 : snapshot%dt = 1.0_dp
173 19078 : IF (PRESENT(dt)) snapshot%dt = dt
174 19078 : IF (wf_history%store_wf) THEN
175 18504 : CALL get_qs_env(qs_env, mos=mos)
176 18504 : IF (.NOT. ASSOCIATED(snapshot%wf)) THEN
177 : CALL fm_pools_create_fm_vect(ao_mo_pools, snapshot%wf, &
178 10440 : name="ws_snap-ws")
179 10440 : CPASSERT(nspins == SIZE(snapshot%wf))
180 : END IF
181 39092 : DO ispin = 1, nspins
182 20588 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
183 39092 : CALL cp_fm_to_fm(mo_coeff, snapshot%wf(ispin))
184 : END DO
185 : ELSE
186 574 : CALL fm_pools_give_back_fm_vect(ao_mo_pools, snapshot%wf)
187 : END IF
188 :
189 19078 : IF (wf_history%store_rho_r) THEN
190 0 : CALL qs_rho_get(rho, rho_r=rho_r)
191 0 : CPASSERT(ASSOCIATED(rho_r))
192 0 : IF (.NOT. ASSOCIATED(snapshot%rho_r)) THEN
193 0 : ALLOCATE (snapshot%rho_r(nspins))
194 0 : DO ispin = 1, nspins
195 0 : CALL auxbas_pw_pool%create_pw(snapshot%rho_r(ispin))
196 : END DO
197 : END IF
198 0 : DO ispin = 1, nspins
199 0 : CALL pw_copy(rho_r(ispin), snapshot%rho_r(ispin))
200 : END DO
201 19078 : ELSE IF (ASSOCIATED(snapshot%rho_r)) THEN
202 0 : DO ispin = 1, SIZE(snapshot%rho_r)
203 0 : CALL auxbas_pw_pool%give_back_pw(snapshot%rho_r(ispin))
204 : END DO
205 0 : DEALLOCATE (snapshot%rho_r)
206 : END IF
207 :
208 19078 : IF (wf_history%store_rho_g) THEN
209 0 : CALL qs_rho_get(rho, rho_g=rho_g)
210 0 : CPASSERT(ASSOCIATED(rho_g))
211 0 : IF (.NOT. ASSOCIATED(snapshot%rho_g)) THEN
212 0 : ALLOCATE (snapshot%rho_g(nspins))
213 0 : DO ispin = 1, nspins
214 0 : CALL auxbas_pw_pool%create_pw(snapshot%rho_g(ispin))
215 : END DO
216 : END IF
217 0 : DO ispin = 1, nspins
218 0 : CALL pw_copy(rho_g(ispin), snapshot%rho_g(ispin))
219 : END DO
220 19078 : ELSE IF (ASSOCIATED(snapshot%rho_g)) THEN
221 0 : DO ispin = 1, SIZE(snapshot%rho_g)
222 0 : CALL auxbas_pw_pool%give_back_pw(snapshot%rho_g(ispin))
223 : END DO
224 0 : DEALLOCATE (snapshot%rho_g)
225 : END IF
226 :
227 19078 : IF (ASSOCIATED(snapshot%rho_ao)) THEN ! the sparsity might be different
228 : ! (future struct:check)
229 270 : CALL dbcsr_deallocate_matrix_set(snapshot%rho_ao)
230 : END IF
231 19078 : IF (wf_history%store_rho_ao) THEN
232 338 : CALL qs_rho_get(rho, rho_ao=rho_ao)
233 338 : CPASSERT(ASSOCIATED(rho_ao))
234 :
235 338 : CALL dbcsr_allocate_matrix_set(snapshot%rho_ao, nspins)
236 836 : DO ispin = 1, nspins
237 498 : ALLOCATE (snapshot%rho_ao(ispin)%matrix)
238 836 : CALL dbcsr_copy(snapshot%rho_ao(ispin)%matrix, rho_ao(ispin)%matrix)
239 : END DO
240 : END IF
241 :
242 19078 : IF (ASSOCIATED(snapshot%rho_ao_kp)) THEN ! the sparsity might be different
243 : ! (future struct:check)
244 214 : CALL dbcsr_deallocate_matrix_set(snapshot%rho_ao_kp)
245 : END IF
246 19078 : IF (wf_history%store_rho_ao_kp) THEN
247 220 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
248 220 : CPASSERT(ASSOCIATED(rho_ao_kp))
249 :
250 220 : nimg = dft_control%nimages
251 220 : CALL dbcsr_allocate_matrix_set(snapshot%rho_ao_kp, nspins, nimg)
252 530 : DO ispin = 1, nspins
253 31010 : DO img = 1, nimg
254 30480 : ALLOCATE (snapshot%rho_ao_kp(ispin, img)%matrix)
255 : CALL dbcsr_copy(snapshot%rho_ao_kp(ispin, img)%matrix, &
256 30790 : rho_ao_kp(ispin, img)%matrix)
257 : END DO
258 : END DO
259 : END IF
260 :
261 19078 : IF (ASSOCIATED(snapshot%overlap)) THEN ! the sparsity might be different
262 : ! (future struct:check)
263 5468 : CALL dbcsr_deallocate_matrix(snapshot%overlap)
264 : END IF
265 19078 : IF (wf_history%store_overlap) THEN
266 14884 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
267 14884 : CPASSERT(ASSOCIATED(matrix_s))
268 14884 : CPASSERT(ASSOCIATED(matrix_s(1)%matrix))
269 14884 : ALLOCATE (snapshot%overlap)
270 14884 : CALL dbcsr_copy(snapshot%overlap, matrix_s(1)%matrix)
271 : END IF
272 :
273 : IF (wf_history%store_frozen_density) THEN
274 : ! do nothing
275 : ! CALL deallocate_matrix_set(snapshot%rho_frozen%rho_ao)
276 : END IF
277 :
278 19078 : CALL timestop(handle)
279 :
280 19078 : END SUBROUTINE wfs_update
281 :
282 : ! **************************************************************************************************
283 : !> \brief ...
284 : !> \param wf_history ...
285 : !> \param interpolation_method_nr the tag of the method used for
286 : !> the extrapolation of the initial density for the next md step
287 : !> (see qs_wf_history_types:wfi_*_method_nr)
288 : !> \param extrapolation_order ...
289 : !> \param has_unit_metric ...
290 : !> \par History
291 : !> 02.2003 created [fawzi]
292 : !> \author fawzi
293 : ! **************************************************************************************************
294 7410 : SUBROUTINE wfi_create(wf_history, interpolation_method_nr, extrapolation_order, &
295 : has_unit_metric)
296 : TYPE(qs_wf_history_type), POINTER :: wf_history
297 : INTEGER, INTENT(in) :: interpolation_method_nr, &
298 : extrapolation_order
299 : LOGICAL, INTENT(IN) :: has_unit_metric
300 :
301 : CHARACTER(len=*), PARAMETER :: routineN = 'wfi_create'
302 :
303 : INTEGER :: handle, i
304 :
305 7410 : CALL timeset(routineN, handle)
306 :
307 7410 : ALLOCATE (wf_history)
308 7410 : wf_history%ref_count = 1
309 7410 : wf_history%memory_depth = 0
310 7410 : wf_history%snapshot_count = 0
311 7410 : wf_history%last_state_index = 1
312 : wf_history%store_wf = .FALSE.
313 : wf_history%store_rho_r = .FALSE.
314 : wf_history%store_rho_g = .FALSE.
315 : wf_history%store_rho_ao = .FALSE.
316 : wf_history%store_rho_ao_kp = .FALSE.
317 : wf_history%store_overlap = .FALSE.
318 : wf_history%store_frozen_density = .FALSE.
319 : NULLIFY (wf_history%past_states)
320 :
321 7410 : wf_history%interpolation_method_nr = interpolation_method_nr
322 :
323 : SELECT CASE (wf_history%interpolation_method_nr)
324 : CASE (wfi_use_guess_method_nr)
325 : wf_history%memory_depth = 0
326 : CASE (wfi_use_prev_wf_method_nr)
327 62 : wf_history%memory_depth = 0
328 : CASE (wfi_use_prev_p_method_nr)
329 62 : wf_history%memory_depth = 1
330 62 : wf_history%store_rho_ao = .TRUE.
331 : CASE (wfi_use_prev_rho_r_method_nr)
332 4 : wf_history%memory_depth = 1
333 4 : wf_history%store_rho_ao = .TRUE.
334 : CASE (wfi_linear_wf_method_nr)
335 4 : wf_history%memory_depth = 2
336 4 : wf_history%store_wf = .TRUE.
337 : CASE (wfi_linear_p_method_nr)
338 4 : wf_history%memory_depth = 2
339 4 : wf_history%store_rho_ao = .TRUE.
340 : CASE (wfi_linear_ps_method_nr)
341 6 : wf_history%memory_depth = 2
342 6 : wf_history%store_wf = .TRUE.
343 6 : IF (.NOT. has_unit_metric) wf_history%store_overlap = .TRUE.
344 : CASE (wfi_ps_method_nr)
345 341 : CALL cite_reference(VandeVondele2005a)
346 341 : wf_history%memory_depth = extrapolation_order + 1
347 341 : wf_history%store_wf = .TRUE.
348 341 : IF (.NOT. has_unit_metric) wf_history%store_overlap = .TRUE.
349 : CASE (wfi_frozen_method_nr)
350 4 : wf_history%memory_depth = 1
351 4 : wf_history%store_frozen_density = .TRUE.
352 : CASE (wfi_aspc_nr)
353 6767 : wf_history%memory_depth = extrapolation_order + 2
354 6767 : wf_history%store_wf = .TRUE.
355 6767 : IF (.NOT. has_unit_metric) wf_history%store_overlap = .TRUE.
356 : CASE default
357 : CALL cp_abort(__LOCATION__, &
358 : "Unknown interpolation method: "// &
359 0 : TRIM(ADJUSTL(cp_to_string(interpolation_method_nr))))
360 7410 : wf_history%interpolation_method_nr = wfi_use_prev_rho_r_method_nr
361 : END SELECT
362 56987 : ALLOCATE (wf_history%past_states(wf_history%memory_depth))
363 :
364 42385 : DO i = 1, SIZE(wf_history%past_states)
365 42385 : NULLIFY (wf_history%past_states(i)%snapshot)
366 : END DO
367 :
368 7410 : CALL timestop(handle)
369 7410 : END SUBROUTINE wfi_create
370 :
371 : ! **************************************************************************************************
372 : !> \brief ...
373 : !> \param wf_history ...
374 : !> \par History
375 : !> 06.2015 created [jhu]
376 : !> \author fawzi
377 : ! **************************************************************************************************
378 170 : SUBROUTINE wfi_create_for_kp(wf_history)
379 : TYPE(qs_wf_history_type), POINTER :: wf_history
380 :
381 170 : CPASSERT(ASSOCIATED(wf_history))
382 170 : IF (wf_history%store_rho_ao) THEN
383 6 : wf_history%store_rho_ao_kp = .TRUE.
384 6 : wf_history%store_rho_ao = .FALSE.
385 : END IF
386 : ! Check for KP compatible methods
387 170 : IF (wf_history%store_wf) THEN
388 0 : CPABORT("WFN based interpolation method not possible for kpoints.")
389 : END IF
390 170 : IF (wf_history%store_frozen_density) THEN
391 0 : CPABORT("Frozen density initialization method not possible for kpoints.")
392 : END IF
393 170 : IF (wf_history%store_overlap) THEN
394 0 : CPABORT("Inconsistent interpolation method for kpoints.")
395 : END IF
396 :
397 170 : END SUBROUTINE wfi_create_for_kp
398 :
399 : ! **************************************************************************************************
400 : !> \brief returns a string describing the interpolation method
401 : !> \param method_nr ...
402 : !> \return ...
403 : !> \par History
404 : !> 02.2003 created [fawzi]
405 : !> \author fawzi
406 : ! **************************************************************************************************
407 10317 : FUNCTION wfi_get_method_label(method_nr) RESULT(res)
408 : INTEGER, INTENT(in) :: method_nr
409 : CHARACTER(len=30) :: res
410 :
411 10317 : res = "unknown"
412 10553 : SELECT CASE (method_nr)
413 : CASE (wfi_use_prev_p_method_nr)
414 236 : res = "previous_p"
415 : CASE (wfi_use_prev_wf_method_nr)
416 210 : res = "previous_wf"
417 : CASE (wfi_use_prev_rho_r_method_nr)
418 4 : res = "previous_rho_r"
419 : CASE (wfi_use_guess_method_nr)
420 3451 : res = "initial_guess"
421 : CASE (wfi_linear_wf_method_nr)
422 2 : res = "mo linear"
423 : CASE (wfi_linear_p_method_nr)
424 2 : res = "P linear"
425 : CASE (wfi_linear_ps_method_nr)
426 6 : res = "PS linear"
427 : CASE (wfi_ps_method_nr)
428 188 : res = "PS Nth order"
429 : CASE (wfi_frozen_method_nr)
430 4 : res = "frozen density approximation"
431 : CASE (wfi_aspc_nr)
432 6214 : res = "ASPC"
433 : CASE default
434 : CALL cp_abort(__LOCATION__, &
435 : "Unknown interpolation method: "// &
436 10317 : TRIM(ADJUSTL(cp_to_string(method_nr))))
437 : END SELECT
438 10317 : END FUNCTION wfi_get_method_label
439 :
440 : ! **************************************************************************************************
441 : !> \brief calculates the new starting state for the scf for the next
442 : !> wf optimization
443 : !> \param wf_history the previous history needed to extrapolate
444 : !> \param qs_env the qs env with the latest result, and that will contain
445 : !> the new starting state
446 : !> \param dt the time at which to extrapolate (wrt. to the last snapshot)
447 : !> \param extrapolation_method_nr returns the extrapolation method used
448 : !> \param orthogonal_wf ...
449 : !> \par History
450 : !> 02.2003 created [fawzi]
451 : !> 11.2003 Joost VandeVondele : Implemented Nth order PS extrapolation
452 : !> \author fawzi
453 : ! **************************************************************************************************
454 20285 : SUBROUTINE wfi_extrapolate(wf_history, qs_env, dt, extrapolation_method_nr, &
455 : orthogonal_wf)
456 : TYPE(qs_wf_history_type), POINTER :: wf_history
457 : TYPE(qs_environment_type), POINTER :: qs_env
458 : REAL(KIND=dp), INTENT(IN) :: dt
459 : INTEGER, INTENT(OUT), OPTIONAL :: extrapolation_method_nr
460 : LOGICAL, INTENT(OUT), OPTIONAL :: orthogonal_wf
461 :
462 : CHARACTER(len=*), PARAMETER :: routineN = 'wfi_extrapolate'
463 :
464 : INTEGER :: actual_extrapolation_method_nr, handle, &
465 : i, img, io_unit, ispin, k, n, nmo, &
466 : nvec, print_level
467 : LOGICAL :: do_kpoints, my_orthogonal_wf, use_overlap
468 : REAL(KIND=dp) :: alpha, t0, t1, t2
469 20285 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools
470 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct, matrix_struct_new
471 : TYPE(cp_fm_type) :: csc, fm_tmp
472 : TYPE(cp_fm_type), POINTER :: mo_coeff
473 : TYPE(cp_logger_type), POINTER :: logger
474 20285 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao, rho_frozen_ao
475 20285 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
476 20285 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
477 : TYPE(qs_rho_type), POINTER :: rho
478 : TYPE(qs_wf_snapshot_type), POINTER :: t0_state, t1_state
479 :
480 20285 : NULLIFY (mos, ao_mo_fm_pools, t0_state, t1_state, mo_coeff, &
481 20285 : rho, rho_ao, rho_frozen_ao)
482 :
483 20285 : use_overlap = wf_history%store_overlap
484 :
485 20285 : CALL timeset(routineN, handle)
486 20285 : logger => cp_get_default_logger()
487 20285 : print_level = logger%iter_info%print_level
488 : io_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%SCF%PRINT%PROGRAM_RUN_INFO", &
489 20285 : extension=".scfLog")
490 :
491 20285 : CPASSERT(ASSOCIATED(wf_history))
492 20285 : CPASSERT(wf_history%ref_count > 0)
493 20285 : CPASSERT(ASSOCIATED(qs_env))
494 20285 : CALL get_qs_env(qs_env, mos=mos, rho=rho, do_kpoints=do_kpoints)
495 20285 : CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
496 : ! chooses the method for this extrapolation
497 20285 : IF (wf_history%snapshot_count < 1) THEN
498 : actual_extrapolation_method_nr = wfi_use_guess_method_nr
499 : ELSE
500 14032 : actual_extrapolation_method_nr = wf_history%interpolation_method_nr
501 : END IF
502 :
503 8 : SELECT CASE (actual_extrapolation_method_nr)
504 : CASE (wfi_linear_wf_method_nr)
505 8 : IF (wf_history%snapshot_count < 2) THEN
506 4 : actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
507 : END IF
508 : CASE (wfi_linear_p_method_nr)
509 8 : IF (wf_history%snapshot_count < 2) THEN
510 4 : IF (do_kpoints) THEN
511 : actual_extrapolation_method_nr = wfi_use_guess_method_nr
512 : ELSE
513 4 : actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
514 : END IF
515 : END IF
516 : CASE (wfi_linear_ps_method_nr)
517 14032 : IF (wf_history%snapshot_count < 2) THEN
518 6 : actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
519 : END IF
520 : END SELECT
521 :
522 20285 : IF (PRESENT(extrapolation_method_nr)) &
523 20285 : extrapolation_method_nr = actual_extrapolation_method_nr
524 20285 : my_orthogonal_wf = .FALSE.
525 :
526 8 : SELECT CASE (actual_extrapolation_method_nr)
527 : CASE (wfi_frozen_method_nr)
528 8 : CPASSERT(.NOT. do_kpoints)
529 8 : t0_state => wfi_get_snapshot(wf_history, wf_index=1)
530 8 : CPASSERT(ASSOCIATED(t0_state%rho_frozen))
531 :
532 8 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
533 8 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
534 :
535 8 : CALL qs_rho_get(t0_state%rho_frozen, rho_ao=rho_frozen_ao)
536 8 : CALL qs_rho_get(rho, rho_ao=rho_ao)
537 16 : DO ispin = 1, SIZE(rho_frozen_ao)
538 : CALL dbcsr_copy(rho_ao(ispin)%matrix, &
539 : rho_frozen_ao(ispin)%matrix, &
540 16 : keep_sparsity=.TRUE.)
541 : END DO
542 : !FM updating rho_ao directly with t0_state%rho_ao would have the
543 : !FM wrong matrix structure
544 8 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
545 8 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
546 :
547 8 : my_orthogonal_wf = .FALSE.
548 : CASE (wfi_use_prev_rho_r_method_nr)
549 8 : t0_state => wfi_get_snapshot(wf_history, wf_index=1)
550 8 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
551 8 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
552 8 : IF (do_kpoints) THEN
553 0 : CPASSERT(ASSOCIATED(t0_state%rho_ao_kp))
554 0 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
555 0 : DO ispin = 1, SIZE(t0_state%rho_ao_kp, 1)
556 0 : DO img = 1, SIZE(t0_state%rho_ao_kp, 2)
557 0 : IF (img > SIZE(rho_ao_kp, 2)) THEN
558 0 : CPWARN("Change in cell neighborlist: might affect quality of initial guess")
559 : ELSE
560 : CALL dbcsr_copy(rho_ao_kp(ispin, img)%matrix, &
561 : t0_state%rho_ao_kp(ispin, img)%matrix, &
562 0 : keep_sparsity=.TRUE.)
563 : END IF
564 : END DO
565 : END DO
566 : ELSE
567 8 : CPASSERT(ASSOCIATED(t0_state%rho_ao))
568 8 : CALL qs_rho_get(rho, rho_ao=rho_ao)
569 16 : DO ispin = 1, SIZE(t0_state%rho_ao)
570 : CALL dbcsr_copy(rho_ao(ispin)%matrix, &
571 : t0_state%rho_ao(ispin)%matrix, &
572 16 : keep_sparsity=.TRUE.)
573 : END DO
574 : END IF
575 : ! Why is rho_g valid at this point ?
576 8 : CALL qs_rho_set(rho, rho_g_valid=.TRUE.)
577 :
578 : ! does nothing
579 : CASE (wfi_use_prev_wf_method_nr)
580 420 : CPASSERT(.NOT. do_kpoints)
581 420 : my_orthogonal_wf = .TRUE.
582 420 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
583 420 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
584 420 : CALL qs_rho_get(rho, rho_ao=rho_ao)
585 888 : DO ispin = 1, SIZE(mos)
586 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, &
587 468 : nmo=nmo)
588 : CALL reorthogonalize_vectors(qs_env, &
589 : v_matrix=mo_coeff, &
590 468 : n_col=nmo)
591 : CALL calculate_density_matrix(mo_set=mos(ispin), &
592 1356 : density_matrix=rho_ao(ispin)%matrix)
593 : END DO
594 420 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
595 420 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
596 :
597 : CASE (wfi_use_prev_p_method_nr)
598 472 : t0_state => wfi_get_snapshot(wf_history, wf_index=1)
599 472 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
600 472 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
601 472 : IF (do_kpoints) THEN
602 214 : CPASSERT(ASSOCIATED(t0_state%rho_ao_kp))
603 214 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
604 516 : DO ispin = 1, SIZE(t0_state%rho_ao_kp, 1)
605 30222 : DO img = 1, SIZE(t0_state%rho_ao_kp, 2)
606 30008 : IF (img > SIZE(rho_ao_kp, 2)) THEN
607 0 : CPWARN("Change in cell neighborlist: might affect quality of initial guess")
608 : ELSE
609 : CALL dbcsr_copy(rho_ao_kp(ispin, img)%matrix, &
610 : t0_state%rho_ao_kp(ispin, img)%matrix, &
611 29706 : keep_sparsity=.TRUE.)
612 : END IF
613 : END DO
614 : END DO
615 : ELSE
616 258 : CPASSERT(ASSOCIATED(t0_state%rho_ao))
617 258 : CALL qs_rho_get(rho, rho_ao=rho_ao)
618 646 : DO ispin = 1, SIZE(t0_state%rho_ao)
619 : CALL dbcsr_copy(rho_ao(ispin)%matrix, &
620 : t0_state%rho_ao(ispin)%matrix, &
621 646 : keep_sparsity=.TRUE.)
622 : END DO
623 : END IF
624 : !FM updating rho_ao directly with t0_state%rho_ao would have the
625 : !FM wrong matrix structure
626 472 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
627 472 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
628 :
629 : CASE (wfi_use_guess_method_nr)
630 : !FM more clean to do it here, but it
631 : !FM might need to read a file (restart) and thus globenv
632 : !FM I do not want globenv here, thus done by the caller
633 : !FM (btw. it also needs the eigensolver, and unless you relocate it
634 : !FM gives circular dependencies)
635 6845 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
636 6845 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
637 : CASE (wfi_linear_wf_method_nr)
638 4 : CPASSERT(.NOT. do_kpoints)
639 4 : t0_state => wfi_get_snapshot(wf_history, wf_index=2)
640 4 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
641 4 : CPASSERT(ASSOCIATED(t0_state))
642 4 : CPASSERT(ASSOCIATED(t1_state))
643 4 : CPASSERT(ASSOCIATED(t0_state%wf))
644 4 : CPASSERT(ASSOCIATED(t1_state%wf))
645 4 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
646 4 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
647 :
648 4 : my_orthogonal_wf = .TRUE.
649 4 : t0 = 0.0_dp
650 4 : t1 = t1_state%dt
651 4 : t2 = t1 + dt
652 4 : CALL qs_rho_get(rho, rho_ao=rho_ao)
653 8 : DO ispin = 1, SIZE(mos)
654 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, &
655 4 : nmo=nmo)
656 : CALL cp_fm_scale_and_add(alpha=0.0_dp, &
657 : matrix_a=mo_coeff, &
658 : matrix_b=t1_state%wf(ispin), &
659 4 : beta=(t2 - t0)/(t1 - t0))
660 : ! this copy should be unnecessary
661 : CALL cp_fm_scale_and_add(alpha=1.0_dp, &
662 : matrix_a=mo_coeff, &
663 4 : beta=(t1 - t2)/(t1 - t0), matrix_b=t0_state%wf(ispin))
664 : CALL reorthogonalize_vectors(qs_env, &
665 : v_matrix=mo_coeff, &
666 4 : n_col=nmo)
667 : CALL calculate_density_matrix(mo_set=mos(ispin), &
668 12 : density_matrix=rho_ao(ispin)%matrix)
669 : END DO
670 4 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
671 :
672 : CALL qs_ks_did_change(qs_env%ks_env, &
673 4 : rho_changed=.TRUE.)
674 : CASE (wfi_linear_p_method_nr)
675 4 : t0_state => wfi_get_snapshot(wf_history, wf_index=2)
676 4 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
677 4 : CPASSERT(ASSOCIATED(t0_state))
678 4 : CPASSERT(ASSOCIATED(t1_state))
679 4 : IF (do_kpoints) THEN
680 0 : CPASSERT(ASSOCIATED(t0_state%rho_ao_kp))
681 0 : CPASSERT(ASSOCIATED(t1_state%rho_ao_kp))
682 : ELSE
683 4 : CPASSERT(ASSOCIATED(t0_state%rho_ao))
684 4 : CPASSERT(ASSOCIATED(t1_state%rho_ao))
685 : END IF
686 4 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
687 4 : CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
688 :
689 4 : t0 = 0.0_dp
690 4 : t1 = t1_state%dt
691 4 : t2 = t1 + dt
692 4 : IF (do_kpoints) THEN
693 0 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
694 0 : DO ispin = 1, SIZE(rho_ao_kp, 1)
695 0 : DO img = 1, SIZE(rho_ao_kp, 2)
696 0 : IF (img > SIZE(t0_state%rho_ao_kp, 2) .OR. &
697 0 : img > SIZE(t1_state%rho_ao_kp, 2)) THEN
698 0 : CPWARN("Change in cell neighborlist: might affect quality of initial guess")
699 : ELSE
700 : CALL dbcsr_add(rho_ao_kp(ispin, img)%matrix, t1_state%rho_ao_kp(ispin, img)%matrix, &
701 0 : alpha_scalar=0.0_dp, beta_scalar=(t2 - t0)/(t1 - t0)) ! this copy should be unnecessary
702 : CALL dbcsr_add(rho_ao_kp(ispin, img)%matrix, t0_state%rho_ao_kp(ispin, img)%matrix, &
703 0 : alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0))
704 : END IF
705 : END DO
706 : END DO
707 : ELSE
708 4 : CALL qs_rho_get(rho, rho_ao=rho_ao)
709 8 : DO ispin = 1, SIZE(rho_ao)
710 : CALL dbcsr_add(rho_ao(ispin)%matrix, t1_state%rho_ao(ispin)%matrix, &
711 4 : alpha_scalar=0.0_dp, beta_scalar=(t2 - t0)/(t1 - t0)) ! this copy should be unnecessary
712 : CALL dbcsr_add(rho_ao(ispin)%matrix, t0_state%rho_ao(ispin)%matrix, &
713 8 : alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0))
714 : END DO
715 : END IF
716 4 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
717 4 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
718 : CASE (wfi_linear_ps_method_nr)
719 : ! wf not calculated, extract with PSC renormalized?
720 : ! use wf_linear?
721 12 : CPASSERT(.NOT. do_kpoints)
722 12 : t0_state => wfi_get_snapshot(wf_history, wf_index=2)
723 12 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
724 12 : CPASSERT(ASSOCIATED(t0_state))
725 12 : CPASSERT(ASSOCIATED(t1_state))
726 12 : CPASSERT(ASSOCIATED(t0_state%wf))
727 12 : CPASSERT(ASSOCIATED(t1_state%wf))
728 12 : IF (wf_history%store_overlap) THEN
729 4 : CPASSERT(ASSOCIATED(t0_state%overlap))
730 4 : CPASSERT(ASSOCIATED(t1_state%overlap))
731 : END IF
732 12 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
733 12 : IF (nvec >= wf_history%memory_depth) THEN
734 12 : IF ((qs_env%scf_control%max_scf_hist .NE. 0) .AND. (qs_env%scf_control%eps_scf_hist .NE. 0)) THEN
735 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
736 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
737 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
738 12 : ELSE IF (qs_env%scf_control%max_scf_hist .NE. 0) THEN
739 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
740 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
741 12 : ELSE IF (qs_env%scf_control%eps_scf_hist .NE. 0) THEN
742 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
743 : END IF
744 : END IF
745 :
746 12 : my_orthogonal_wf = .TRUE.
747 : ! use PS_2=2 PS_1-PS_0
748 : ! C_2 comes from using PS_2 as a projector acting on C_1
749 12 : CALL qs_rho_get(rho, rho_ao=rho_ao)
750 24 : DO ispin = 1, SIZE(mos)
751 12 : NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
752 12 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
753 : CALL cp_fm_get_info(mo_coeff, nrow_global=n, ncol_global=k, &
754 12 : matrix_struct=matrix_struct)
755 : CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, &
756 12 : nrow_global=k, ncol_global=k)
757 12 : CALL cp_fm_create(csc, matrix_struct_new)
758 12 : CALL cp_fm_struct_release(matrix_struct_new)
759 :
760 12 : IF (use_overlap) THEN
761 4 : CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), mo_coeff, k)
762 4 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), mo_coeff, 0.0_dp, csc)
763 : ELSE
764 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
765 8 : t1_state%wf(ispin), 0.0_dp, csc)
766 : END IF
767 12 : CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, mo_coeff)
768 12 : CALL cp_fm_release(csc)
769 12 : CALL cp_fm_scale_and_add(-1.0_dp, mo_coeff, 2.0_dp, t1_state%wf(ispin))
770 : CALL reorthogonalize_vectors(qs_env, &
771 : v_matrix=mo_coeff, &
772 12 : n_col=k)
773 : CALL calculate_density_matrix(mo_set=mos(ispin), &
774 48 : density_matrix=rho_ao(ispin)%matrix)
775 : END DO
776 12 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
777 12 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
778 :
779 : CASE (wfi_ps_method_nr)
780 376 : CPASSERT(.NOT. do_kpoints)
781 : ! figure out the actual number of vectors to use in the extrapolation:
782 376 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
783 376 : CPASSERT(nvec .GT. 0)
784 376 : IF (nvec >= wf_history%memory_depth) THEN
785 178 : IF ((qs_env%scf_control%max_scf_hist .NE. 0) .AND. (qs_env%scf_control%eps_scf_hist .NE. 0)) THEN
786 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
787 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
788 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
789 178 : ELSE IF (qs_env%scf_control%max_scf_hist .NE. 0) THEN
790 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
791 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
792 178 : ELSE IF (qs_env%scf_control%eps_scf_hist .NE. 0) THEN
793 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
794 : END IF
795 : END IF
796 :
797 376 : my_orthogonal_wf = .TRUE.
798 830 : DO ispin = 1, SIZE(mos)
799 454 : NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
800 454 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
801 : CALL cp_fm_get_info(mo_coeff, nrow_global=n, ncol_global=k, &
802 454 : matrix_struct=matrix_struct)
803 454 : CALL cp_fm_create(fm_tmp, matrix_struct)
804 : CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, &
805 454 : nrow_global=k, ncol_global=k)
806 454 : CALL cp_fm_create(csc, matrix_struct_new)
807 454 : CALL cp_fm_struct_release(matrix_struct_new)
808 : ! first the most recent
809 454 : t1_state => wfi_get_snapshot(wf_history, wf_index=1)
810 454 : CALL cp_fm_to_fm(t1_state%wf(ispin), mo_coeff)
811 454 : alpha = nvec
812 454 : CALL cp_fm_scale(alpha, mo_coeff)
813 454 : CALL qs_rho_get(rho, rho_ao=rho_ao)
814 968 : DO i = 2, nvec
815 514 : t0_state => wfi_get_snapshot(wf_history, wf_index=i)
816 514 : IF (use_overlap) THEN
817 476 : CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
818 476 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
819 : ELSE
820 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
821 38 : t1_state%wf(ispin), 0.0_dp, csc)
822 : END IF
823 514 : CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
824 514 : alpha = -1.0_dp*alpha*REAL(nvec - i + 1, dp)/REAL(i, dp)
825 968 : CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp)
826 : END DO
827 :
828 454 : CALL cp_fm_release(csc)
829 454 : CALL cp_fm_release(fm_tmp)
830 : CALL reorthogonalize_vectors(qs_env, &
831 : v_matrix=mo_coeff, &
832 454 : n_col=k)
833 : CALL calculate_density_matrix(mo_set=mos(ispin), &
834 1738 : density_matrix=rho_ao(ispin)%matrix)
835 : END DO
836 376 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
837 376 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
838 :
839 : CASE (wfi_aspc_nr)
840 12136 : CPASSERT(.NOT. do_kpoints)
841 12136 : CALL cite_reference(Kolafa2004)
842 12136 : CALL cite_reference(Kuhne2007)
843 : ! figure out the actual number of vectors to use in the extrapolation:
844 12136 : nvec = MIN(wf_history%memory_depth, wf_history%snapshot_count)
845 12136 : CPASSERT(nvec .GT. 0)
846 12136 : IF (nvec >= wf_history%memory_depth) THEN
847 7870 : IF ((qs_env%scf_control%max_scf_hist .NE. 0) .AND. &
848 : (qs_env%scf_control%eps_scf_hist .NE. 0)) THEN
849 18 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
850 18 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
851 18 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
852 7852 : ELSE IF (qs_env%scf_control%max_scf_hist .NE. 0) THEN
853 62 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
854 62 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
855 7790 : ELSE IF (qs_env%scf_control%eps_scf_hist .NE. 0) THEN
856 8 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
857 : END IF
858 : END IF
859 :
860 12136 : my_orthogonal_wf = .TRUE.
861 12136 : CALL qs_rho_get(rho, rho_ao=rho_ao)
862 24833 : DO ispin = 1, SIZE(mos)
863 12697 : NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
864 12697 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
865 : CALL cp_fm_get_info(mo_coeff, &
866 : nrow_global=n, &
867 : ncol_global=k, &
868 12697 : matrix_struct=matrix_struct)
869 12697 : CALL cp_fm_create(fm_tmp, matrix_struct)
870 : CALL cp_fm_struct_create(matrix_struct_new, &
871 : template_fmstruct=matrix_struct, &
872 : nrow_global=k, &
873 12697 : ncol_global=k)
874 12697 : CALL cp_fm_create(csc, matrix_struct_new)
875 12697 : CALL cp_fm_struct_release(matrix_struct_new)
876 : ! first the most recent
877 : t1_state => wfi_get_snapshot(wf_history, &
878 12697 : wf_index=1)
879 12697 : CALL cp_fm_to_fm(t1_state%wf(ispin), mo_coeff)
880 12697 : alpha = REAL(4*nvec - 2, KIND=dp)/REAL(nvec + 1, KIND=dp)
881 12697 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
882 : WRITE (UNIT=io_unit, FMT="(/,T2,A,/,/,T3,A,I0,/,/,T3,A2,I0,A4,F10.6)") &
883 2360 : "Parameters for the always stable predictor-corrector (ASPC) method:", &
884 2360 : "ASPC order: ", MAX(nvec - 2, 0), &
885 4720 : "B(", 1, ") = ", alpha
886 : END IF
887 12697 : CALL cp_fm_scale(alpha, mo_coeff)
888 :
889 49367 : DO i = 2, nvec
890 36670 : t0_state => wfi_get_snapshot(wf_history, wf_index=i)
891 36670 : IF (use_overlap) THEN
892 25190 : CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
893 25190 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), fm_tmp, 0.0_dp, csc)
894 : ELSE
895 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, t0_state%wf(ispin), &
896 11480 : t1_state%wf(ispin), 0.0_dp, csc)
897 : END IF
898 36670 : CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
899 : alpha = (-1.0_dp)**(i + 1)*REAL(i, KIND=dp)* &
900 36670 : binomial(2*nvec, nvec - i)/binomial(2*nvec - 2, nvec - 1)
901 36670 : IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
902 : WRITE (UNIT=io_unit, FMT="(T3,A2,I0,A4,F10.6)") &
903 6843 : "B(", i, ") = ", alpha
904 : END IF
905 49367 : CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp)
906 : END DO
907 12697 : CALL cp_fm_release(csc)
908 12697 : CALL cp_fm_release(fm_tmp)
909 : CALL reorthogonalize_vectors(qs_env, &
910 : v_matrix=mo_coeff, &
911 12697 : n_col=k)
912 : CALL calculate_density_matrix(mo_set=mos(ispin), &
913 50227 : density_matrix=rho_ao(ispin)%matrix)
914 : END DO
915 12136 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
916 12136 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
917 :
918 : CASE default
919 : CALL cp_abort(__LOCATION__, &
920 : "Unknown interpolation method: "// &
921 20285 : TRIM(ADJUSTL(cp_to_string(wf_history%interpolation_method_nr))))
922 : END SELECT
923 20285 : IF (PRESENT(orthogonal_wf)) orthogonal_wf = my_orthogonal_wf
924 : CALL cp_print_key_finished_output(io_unit, logger, qs_env%input, &
925 20285 : "DFT%SCF%PRINT%PROGRAM_RUN_INFO")
926 20285 : CALL timestop(handle)
927 20285 : END SUBROUTINE wfi_extrapolate
928 :
929 : ! **************************************************************************************************
930 : !> \brief Decides if scf control variables has to changed due
931 : !> to using a WF extrapolation.
932 : !> \param qs_env The QS environment
933 : !> \param nvec ...
934 : !> \par History
935 : !> 11.2006 created [TdK]
936 : !> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch)
937 : ! **************************************************************************************************
938 7761 : ELEMENTAL SUBROUTINE wfi_set_history_variables(qs_env, nvec)
939 : TYPE(qs_environment_type), INTENT(INOUT) :: qs_env
940 : INTEGER, INTENT(IN) :: nvec
941 :
942 7761 : IF (nvec >= qs_env%wf_history%memory_depth) THEN
943 1709 : IF ((qs_env%scf_control%max_scf_hist .NE. 0) .AND. (qs_env%scf_control%eps_scf_hist .NE. 0)) THEN
944 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
945 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
946 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
947 1709 : ELSE IF (qs_env%scf_control%max_scf_hist .NE. 0) THEN
948 0 : qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
949 0 : qs_env%scf_control%outer_scf%have_scf = .FALSE.
950 1709 : ELSE IF (qs_env%scf_control%eps_scf_hist .NE. 0) THEN
951 0 : qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
952 0 : qs_env%scf_control%outer_scf%eps_scf = qs_env%scf_control%eps_scf_hist
953 : END IF
954 : END IF
955 :
956 7761 : END SUBROUTINE wfi_set_history_variables
957 :
958 : ! **************************************************************************************************
959 : !> \brief updates the snapshot buffer, taking a new snapshot
960 : !> \param wf_history the history buffer to update
961 : !> \param qs_env the qs_env we get the info from
962 : !> \param dt ...
963 : !> \par History
964 : !> 02.2003 created [fawzi]
965 : !> \author fawzi
966 : ! **************************************************************************************************
967 20289 : SUBROUTINE wfi_update(wf_history, qs_env, dt)
968 : TYPE(qs_wf_history_type), POINTER :: wf_history
969 : TYPE(qs_environment_type), POINTER :: qs_env
970 : REAL(KIND=dp), INTENT(in) :: dt
971 :
972 20289 : CPASSERT(ASSOCIATED(wf_history))
973 20289 : CPASSERT(wf_history%ref_count > 0)
974 20289 : CPASSERT(ASSOCIATED(qs_env))
975 :
976 20289 : wf_history%snapshot_count = wf_history%snapshot_count + 1
977 20289 : IF (wf_history%memory_depth > 0) THEN
978 : wf_history%last_state_index = MODULO(wf_history%snapshot_count, &
979 19078 : wf_history%memory_depth) + 1
980 : CALL wfs_update(snapshot=wf_history%past_states &
981 : (wf_history%last_state_index)%snapshot, wf_history=wf_history, &
982 19078 : qs_env=qs_env, dt=dt)
983 : END IF
984 20289 : END SUBROUTINE wfi_update
985 :
986 : ! **************************************************************************************************
987 : !> \brief reorthogonalizes the mos
988 : !> \param qs_env the qs_env in which to orthogonalize
989 : !> \param v_matrix the vectors to orthogonalize
990 : !> \param n_col number of column of v to orthogonalize
991 : !> \par History
992 : !> 04.2003 created [fawzi]
993 : !> \author Fawzi Mohamed
994 : ! **************************************************************************************************
995 27270 : SUBROUTINE reorthogonalize_vectors(qs_env, v_matrix, n_col)
996 : TYPE(qs_environment_type), POINTER :: qs_env
997 : TYPE(cp_fm_type), INTENT(IN) :: v_matrix
998 : INTEGER, INTENT(in), OPTIONAL :: n_col
999 :
1000 : CHARACTER(len=*), PARAMETER :: routineN = 'reorthogonalize_vectors'
1001 :
1002 : INTEGER :: handle, my_n_col
1003 : LOGICAL :: has_unit_metric, &
1004 : ortho_contains_cholesky, &
1005 : smearing_is_used
1006 : TYPE(cp_fm_pool_type), POINTER :: maxao_maxmo_fm_pool
1007 13635 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1008 : TYPE(dft_control_type), POINTER :: dft_control
1009 : TYPE(qs_matrix_pools_type), POINTER :: mpools
1010 : TYPE(qs_scf_env_type), POINTER :: scf_env
1011 : TYPE(scf_control_type), POINTER :: scf_control
1012 :
1013 13635 : NULLIFY (scf_env, scf_control, maxao_maxmo_fm_pool, matrix_s, mpools, dft_control)
1014 13635 : CALL timeset(routineN, handle)
1015 :
1016 13635 : CPASSERT(ASSOCIATED(qs_env))
1017 :
1018 13635 : CALL cp_fm_get_info(v_matrix, ncol_global=my_n_col)
1019 13635 : IF (PRESENT(n_col)) my_n_col = n_col
1020 : CALL get_qs_env(qs_env, mpools=mpools, &
1021 : scf_env=scf_env, &
1022 : scf_control=scf_control, &
1023 : matrix_s=matrix_s, &
1024 13635 : dft_control=dft_control)
1025 13635 : CALL mpools_get(mpools, maxao_maxmo_fm_pool=maxao_maxmo_fm_pool)
1026 13635 : IF (ASSOCIATED(scf_env)) THEN
1027 : ortho_contains_cholesky = (scf_env%method /= ot_method_nr) .AND. &
1028 : (scf_env%cholesky_method > 0) .AND. &
1029 13635 : ASSOCIATED(scf_env%ortho)
1030 : ELSE
1031 : ortho_contains_cholesky = .FALSE.
1032 : END IF
1033 :
1034 13635 : CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
1035 13635 : smearing_is_used = .FALSE.
1036 13635 : IF (dft_control%smear) THEN
1037 514 : smearing_is_used = .TRUE.
1038 : END IF
1039 :
1040 13635 : IF (has_unit_metric) THEN
1041 3390 : CALL make_basis_simple(v_matrix, my_n_col)
1042 10245 : ELSE IF (smearing_is_used) THEN
1043 : CALL make_basis_lowdin(vmatrix=v_matrix, ncol=my_n_col, &
1044 514 : matrix_s=matrix_s(1)%matrix)
1045 9731 : ELSE IF (ortho_contains_cholesky) THEN
1046 : CALL make_basis_cholesky(vmatrix=v_matrix, ncol=my_n_col, &
1047 6606 : ortho=scf_env%ortho)
1048 : ELSE
1049 3125 : CALL make_basis_sm(v_matrix, my_n_col, matrix_s(1)%matrix)
1050 : END IF
1051 13635 : CALL timestop(handle)
1052 13635 : END SUBROUTINE reorthogonalize_vectors
1053 :
1054 : ! **************************************************************************************************
1055 : !> \brief purges wf_history retaining only the latest snapshot
1056 : !> \param qs_env the qs env with the latest result, and that will contain
1057 : !> the purged wf_history
1058 : !> \par History
1059 : !> 05.2016 created [Nico Holmberg]
1060 : !> \author Nico Holmberg
1061 : ! **************************************************************************************************
1062 0 : SUBROUTINE wfi_purge_history(qs_env)
1063 : TYPE(qs_environment_type), POINTER :: qs_env
1064 :
1065 : CHARACTER(len=*), PARAMETER :: routineN = 'wfi_purge_history'
1066 :
1067 : INTEGER :: handle, io_unit, print_level
1068 : TYPE(cp_logger_type), POINTER :: logger
1069 : TYPE(dft_control_type), POINTER :: dft_control
1070 : TYPE(qs_wf_history_type), POINTER :: wf_history
1071 :
1072 0 : NULLIFY (dft_control, wf_history)
1073 :
1074 0 : CALL timeset(routineN, handle)
1075 0 : logger => cp_get_default_logger()
1076 0 : print_level = logger%iter_info%print_level
1077 : io_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%SCF%PRINT%PROGRAM_RUN_INFO", &
1078 0 : extension=".scfLog")
1079 :
1080 0 : CPASSERT(ASSOCIATED(qs_env))
1081 0 : CPASSERT(ASSOCIATED(qs_env%wf_history))
1082 0 : CPASSERT(qs_env%wf_history%ref_count > 0)
1083 0 : CALL get_qs_env(qs_env, dft_control=dft_control)
1084 :
1085 0 : SELECT CASE (qs_env%wf_history%interpolation_method_nr)
1086 : CASE (wfi_use_guess_method_nr, wfi_use_prev_wf_method_nr, &
1087 : wfi_use_prev_p_method_nr, wfi_use_prev_rho_r_method_nr, &
1088 : wfi_frozen_method_nr)
1089 : ! do nothing
1090 : CASE (wfi_linear_wf_method_nr, wfi_linear_p_method_nr, &
1091 : wfi_linear_ps_method_nr, wfi_ps_method_nr, &
1092 : wfi_aspc_nr)
1093 0 : IF (qs_env%wf_history%snapshot_count .GE. 2) THEN
1094 0 : IF (debug_this_module .AND. io_unit > 0) &
1095 0 : WRITE (io_unit, FMT="(T2,A)") "QS| Purging WFN history"
1096 : CALL wfi_create(wf_history, interpolation_method_nr= &
1097 : dft_control%qs_control%wf_interpolation_method_nr, &
1098 : extrapolation_order=dft_control%qs_control%wf_extrapolation_order, &
1099 0 : has_unit_metric=qs_env%has_unit_metric)
1100 : CALL set_qs_env(qs_env=qs_env, &
1101 0 : wf_history=wf_history)
1102 0 : CALL wfi_release(wf_history)
1103 0 : CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
1104 : END IF
1105 : CASE DEFAULT
1106 0 : CPABORT("Unknown extrapolation method.")
1107 : END SELECT
1108 0 : CALL timestop(handle)
1109 :
1110 0 : END SUBROUTINE wfi_purge_history
1111 :
1112 : END MODULE qs_wf_history_methods
|
