(git:ccc2433)
qs_wf_history_methods.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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 ! **************************************************************************************************
22  USE bibliography, ONLY: kolafa2004,&
23  kuhne2007,&
25  cite_reference
26  USE cp_control_types, ONLY: dft_control_type
30  USE cp_fm_basic_linalg, ONLY: cp_fm_scale,&
32  USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
33  cp_fm_pool_type,&
34  fm_pools_create_fm_vect,&
35  fm_pools_give_back_fm_vect
38  cp_fm_struct_type
39  USE cp_fm_types, ONLY: cp_fm_create,&
41  cp_fm_release,&
42  cp_fm_to_fm,&
43  cp_fm_type
45  cp_logger_type,&
46  cp_to_string
50  USE dbcsr_api, ONLY: dbcsr_add,&
51  dbcsr_copy,&
52  dbcsr_deallocate_matrix,&
53  dbcsr_p_type
54  USE input_constants, ONLY: &
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,&
71  USE qs_ks_types, ONLY: qs_ks_did_change
72  USE qs_matrix_pools, ONLY: mpools_get,&
73  qs_matrix_pools_type
78  USE qs_mo_types, ONLY: get_mo_set,&
79  mo_set_type
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,&
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 
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  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  snapshot%dt = 1.0_dp
120  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  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  TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_pools
143  TYPE(cp_fm_type), POINTER :: mo_coeff
144  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
145  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
146  TYPE(dft_control_type), POINTER :: dft_control
147  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
148  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  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
152  TYPE(qs_rho_type), POINTER :: rho
153 
154  CALL timeset(routinen, handle)
155 
156  NULLIFY (pw_env, auxbas_pw_pool, ao_mo_pools, dft_control, mos, mo_coeff, &
157  rho, rho_r, rho_g, rho_ao, matrix_s)
158  CALL get_qs_env(qs_env, pw_env=pw_env, &
159  dft_control=dft_control, rho=rho)
160  CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_pools)
161  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
162 
163  cpassert(ASSOCIATED(wf_history))
164  cpassert(ASSOCIATED(dft_control))
165  IF (.NOT. ASSOCIATED(snapshot)) THEN
166  ALLOCATE (snapshot)
167  CALL wfs_create(snapshot)
168  END IF
169  cpassert(wf_history%ref_count > 0)
170 
171  nspins = dft_control%nspins
172  snapshot%dt = 1.0_dp
173  IF (PRESENT(dt)) snapshot%dt = dt
174  IF (wf_history%store_wf) THEN
175  CALL get_qs_env(qs_env, mos=mos)
176  IF (.NOT. ASSOCIATED(snapshot%wf)) THEN
177  CALL fm_pools_create_fm_vect(ao_mo_pools, snapshot%wf, &
178  name="ws_snap-ws")
179  cpassert(nspins == SIZE(snapshot%wf))
180  END IF
181  DO ispin = 1, nspins
182  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
183  CALL cp_fm_to_fm(mo_coeff, snapshot%wf(ispin))
184  END DO
185  ELSE
186  CALL fm_pools_give_back_fm_vect(ao_mo_pools, snapshot%wf)
187  END IF
188 
189  IF (wf_history%store_rho_r) THEN
190  CALL qs_rho_get(rho, rho_r=rho_r)
191  cpassert(ASSOCIATED(rho_r))
192  IF (.NOT. ASSOCIATED(snapshot%rho_r)) THEN
193  ALLOCATE (snapshot%rho_r(nspins))
194  DO ispin = 1, nspins
195  CALL auxbas_pw_pool%create_pw(snapshot%rho_r(ispin))
196  END DO
197  END IF
198  DO ispin = 1, nspins
199  CALL pw_copy(rho_r(ispin), snapshot%rho_r(ispin))
200  END DO
201  ELSE IF (ASSOCIATED(snapshot%rho_r)) THEN
202  DO ispin = 1, SIZE(snapshot%rho_r)
203  CALL auxbas_pw_pool%give_back_pw(snapshot%rho_r(ispin))
204  END DO
205  DEALLOCATE (snapshot%rho_r)
206  END IF
207 
208  IF (wf_history%store_rho_g) THEN
209  CALL qs_rho_get(rho, rho_g=rho_g)
210  cpassert(ASSOCIATED(rho_g))
211  IF (.NOT. ASSOCIATED(snapshot%rho_g)) THEN
212  ALLOCATE (snapshot%rho_g(nspins))
213  DO ispin = 1, nspins
214  CALL auxbas_pw_pool%create_pw(snapshot%rho_g(ispin))
215  END DO
216  END IF
217  DO ispin = 1, nspins
218  CALL pw_copy(rho_g(ispin), snapshot%rho_g(ispin))
219  END DO
220  ELSE IF (ASSOCIATED(snapshot%rho_g)) THEN
221  DO ispin = 1, SIZE(snapshot%rho_g)
222  CALL auxbas_pw_pool%give_back_pw(snapshot%rho_g(ispin))
223  END DO
224  DEALLOCATE (snapshot%rho_g)
225  END IF
226 
227  IF (ASSOCIATED(snapshot%rho_ao)) THEN ! the sparsity might be different
228  ! (future struct:check)
229  CALL dbcsr_deallocate_matrix_set(snapshot%rho_ao)
230  END IF
231  IF (wf_history%store_rho_ao) THEN
232  CALL qs_rho_get(rho, rho_ao=rho_ao)
233  cpassert(ASSOCIATED(rho_ao))
234 
235  CALL dbcsr_allocate_matrix_set(snapshot%rho_ao, nspins)
236  DO ispin = 1, nspins
237  ALLOCATE (snapshot%rho_ao(ispin)%matrix)
238  CALL dbcsr_copy(snapshot%rho_ao(ispin)%matrix, rho_ao(ispin)%matrix)
239  END DO
240  END IF
241 
242  IF (ASSOCIATED(snapshot%rho_ao_kp)) THEN ! the sparsity might be different
243  ! (future struct:check)
244  CALL dbcsr_deallocate_matrix_set(snapshot%rho_ao_kp)
245  END IF
246  IF (wf_history%store_rho_ao_kp) THEN
247  CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
248  cpassert(ASSOCIATED(rho_ao_kp))
249 
250  nimg = dft_control%nimages
251  CALL dbcsr_allocate_matrix_set(snapshot%rho_ao_kp, nspins, nimg)
252  DO ispin = 1, nspins
253  DO img = 1, nimg
254  ALLOCATE (snapshot%rho_ao_kp(ispin, img)%matrix)
255  CALL dbcsr_copy(snapshot%rho_ao_kp(ispin, img)%matrix, &
256  rho_ao_kp(ispin, img)%matrix)
257  END DO
258  END DO
259  END IF
260 
261  IF (ASSOCIATED(snapshot%overlap)) THEN ! the sparsity might be different
262  ! (future struct:check)
263  CALL dbcsr_deallocate_matrix(snapshot%overlap)
264  END IF
265  IF (wf_history%store_overlap) THEN
266  CALL get_qs_env(qs_env, matrix_s=matrix_s)
267  cpassert(ASSOCIATED(matrix_s))
268  cpassert(ASSOCIATED(matrix_s(1)%matrix))
269  ALLOCATE (snapshot%overlap)
270  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  CALL timestop(handle)
279 
280  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  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  CALL timeset(routinen, handle)
306 
307  ALLOCATE (wf_history)
308  wf_history%ref_count = 1
309  wf_history%memory_depth = 0
310  wf_history%snapshot_count = 0
311  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  wf_history%interpolation_method_nr = interpolation_method_nr
322 
323  SELECT CASE (wf_history%interpolation_method_nr)
325  wf_history%memory_depth = 0
327  wf_history%memory_depth = 0
329  wf_history%memory_depth = 1
330  wf_history%store_rho_ao = .true.
332  wf_history%memory_depth = 1
333  wf_history%store_rho_ao = .true.
335  wf_history%memory_depth = 2
336  wf_history%store_wf = .true.
338  wf_history%memory_depth = 2
339  wf_history%store_rho_ao = .true.
341  wf_history%memory_depth = 2
342  wf_history%store_wf = .true.
343  IF (.NOT. has_unit_metric) wf_history%store_overlap = .true.
344  CASE (wfi_ps_method_nr)
345  CALL cite_reference(vandevondele2005a)
346  wf_history%memory_depth = extrapolation_order + 1
347  wf_history%store_wf = .true.
348  IF (.NOT. has_unit_metric) wf_history%store_overlap = .true.
349  CASE (wfi_frozen_method_nr)
350  wf_history%memory_depth = 1
351  wf_history%store_frozen_density = .true.
352  CASE (wfi_aspc_nr)
353  wf_history%memory_depth = extrapolation_order + 2
354  wf_history%store_wf = .true.
355  IF (.NOT. has_unit_metric) wf_history%store_overlap = .true.
356  CASE default
357  CALL cp_abort(__location__, &
358  "Unknown interpolation method: "// &
359  trim(adjustl(cp_to_string(interpolation_method_nr))))
360  wf_history%interpolation_method_nr = wfi_use_prev_rho_r_method_nr
361  END SELECT
362  ALLOCATE (wf_history%past_states(wf_history%memory_depth))
363 
364  DO i = 1, SIZE(wf_history%past_states)
365  NULLIFY (wf_history%past_states(i)%snapshot)
366  END DO
367 
368  CALL timestop(handle)
369  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  SUBROUTINE wfi_create_for_kp(wf_history)
379  TYPE(qs_wf_history_type), POINTER :: wf_history
380 
381  cpassert(ASSOCIATED(wf_history))
382  IF (wf_history%store_rho_ao) THEN
383  wf_history%store_rho_ao_kp = .true.
384  wf_history%store_rho_ao = .false.
385  END IF
386  ! Check for KP compatible methods
387  IF (wf_history%store_wf) THEN
388  cpabort("WFN based interpolation method not possible for kpoints.")
389  END IF
390  IF (wf_history%store_frozen_density) THEN
391  cpabort("Frozen density initialization method not possible for kpoints.")
392  END IF
393  IF (wf_history%store_overlap) THEN
394  cpabort("Inconsistent interpolation method for kpoints.")
395  END IF
396 
397  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  FUNCTION wfi_get_method_label(method_nr) RESULT(res)
408  INTEGER, INTENT(in) :: method_nr
409  CHARACTER(len=30) :: res
410 
411  res = "unknown"
412  SELECT CASE (method_nr)
414  res = "previous_p"
416  res = "previous_wf"
418  res = "previous_rho_r"
420  res = "initial_guess"
422  res = "mo linear"
424  res = "P linear"
426  res = "PS linear"
427  CASE (wfi_ps_method_nr)
428  res = "PS Nth order"
429  CASE (wfi_frozen_method_nr)
430  res = "frozen density approximation"
431  CASE (wfi_aspc_nr)
432  res = "ASPC"
433  CASE default
434  CALL cp_abort(__location__, &
435  "Unknown interpolation method: "// &
436  trim(adjustl(cp_to_string(method_nr))))
437  END SELECT
438  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  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  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  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao, rho_frozen_ao
475  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
476  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  NULLIFY (mos, ao_mo_fm_pools, t0_state, t1_state, mo_coeff, &
481  rho, rho_ao, rho_frozen_ao)
482 
483  use_overlap = wf_history%store_overlap
484 
485  CALL timeset(routinen, handle)
486  logger => cp_get_default_logger()
487  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  extension=".scfLog")
490 
491  cpassert(ASSOCIATED(wf_history))
492  cpassert(wf_history%ref_count > 0)
493  cpassert(ASSOCIATED(qs_env))
494  CALL get_qs_env(qs_env, mos=mos, rho=rho, do_kpoints=do_kpoints)
495  CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
496  ! chooses the method for this extrapolation
497  IF (wf_history%snapshot_count < 1) THEN
498  actual_extrapolation_method_nr = wfi_use_guess_method_nr
499  ELSE
500  actual_extrapolation_method_nr = wf_history%interpolation_method_nr
501  END IF
502 
503  SELECT CASE (actual_extrapolation_method_nr)
505  IF (wf_history%snapshot_count < 2) THEN
506  actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
507  END IF
509  IF (wf_history%snapshot_count < 2) THEN
510  IF (do_kpoints) THEN
511  actual_extrapolation_method_nr = wfi_use_guess_method_nr
512  ELSE
513  actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
514  END IF
515  END IF
517  IF (wf_history%snapshot_count < 2) THEN
518  actual_extrapolation_method_nr = wfi_use_prev_wf_method_nr
519  END IF
520  END SELECT
521 
522  IF (PRESENT(extrapolation_method_nr)) &
523  extrapolation_method_nr = actual_extrapolation_method_nr
524  my_orthogonal_wf = .false.
525 
526  SELECT CASE (actual_extrapolation_method_nr)
527  CASE (wfi_frozen_method_nr)
528  cpassert(.NOT. do_kpoints)
529  t0_state => wfi_get_snapshot(wf_history, wf_index=1)
530  cpassert(ASSOCIATED(t0_state%rho_frozen))
531 
532  nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
533  CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
534 
535  CALL qs_rho_get(t0_state%rho_frozen, rho_ao=rho_frozen_ao)
536  CALL qs_rho_get(rho, rho_ao=rho_ao)
537  DO ispin = 1, SIZE(rho_frozen_ao)
538  CALL dbcsr_copy(rho_ao(ispin)%matrix, &
539  rho_frozen_ao(ispin)%matrix, &
540  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  CALL qs_rho_update_rho(rho, qs_env=qs_env)
545  CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
546 
547  my_orthogonal_wf = .false.
549  t0_state => wfi_get_snapshot(wf_history, wf_index=1)
550  nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
551  CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
552  IF (do_kpoints) THEN
553  cpassert(ASSOCIATED(t0_state%rho_ao_kp))
554  CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
555  DO ispin = 1, SIZE(t0_state%rho_ao_kp, 1)
556  DO img = 1, SIZE(t0_state%rho_ao_kp, 2)
557  IF (img > SIZE(rho_ao_kp, 2)) THEN
558  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  keep_sparsity=.true.)
563  END IF
564  END DO
565  END DO
566  ELSE
567  cpassert(ASSOCIATED(t0_state%rho_ao))
568  CALL qs_rho_get(rho, rho_ao=rho_ao)
569  DO ispin = 1, SIZE(t0_state%rho_ao)
570  CALL dbcsr_copy(rho_ao(ispin)%matrix, &
571  t0_state%rho_ao(ispin)%matrix, &
572  keep_sparsity=.true.)
573  END DO
574  END IF
575  ! Why is rho_g valid at this point ?
576  CALL qs_rho_set(rho, rho_g_valid=.true.)
577 
578  ! does nothing
580  cpassert(.NOT. do_kpoints)
581  my_orthogonal_wf = .true.
582  nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
583  CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
584  CALL qs_rho_get(rho, rho_ao=rho_ao)
585  DO ispin = 1, SIZE(mos)
586  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, &
587  nmo=nmo)
588  CALL reorthogonalize_vectors(qs_env, &
589  v_matrix=mo_coeff, &
590  n_col=nmo)
591  CALL calculate_density_matrix(mo_set=mos(ispin), &
592  density_matrix=rho_ao(ispin)%matrix)
593  END DO
594  CALL qs_rho_update_rho(rho, qs_env=qs_env)
595  CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
596 
598  t0_state => wfi_get_snapshot(wf_history, wf_index=1)
599  nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
600  CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
601  IF (do_kpoints) THEN
602  cpassert(ASSOCIATED(t0_state%rho_ao_kp))
603  CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
604  DO ispin = 1, SIZE(t0_state%rho_ao_kp, 1)
605  DO img = 1, SIZE(t0_state%rho_ao_kp, 2)
606  IF (img > SIZE(rho_ao_kp, 2)) THEN
607  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  keep_sparsity=.true.)
612  END IF
613  END DO
614  END DO
615  ELSE
616  cpassert(ASSOCIATED(t0_state%rho_ao))
617  CALL qs_rho_get(rho, rho_ao=rho_ao)
618  DO ispin = 1, SIZE(t0_state%rho_ao)
619  CALL dbcsr_copy(rho_ao(ispin)%matrix, &
620  t0_state%rho_ao(ispin)%matrix, &
621  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  CALL qs_rho_update_rho(rho, qs_env=qs_env)
627  CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
628 
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  nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
636  CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
638  cpassert(.NOT. do_kpoints)
639  t0_state => wfi_get_snapshot(wf_history, wf_index=2)
640  t1_state => wfi_get_snapshot(wf_history, wf_index=1)
641  cpassert(ASSOCIATED(t0_state))
642  cpassert(ASSOCIATED(t1_state))
643  cpassert(ASSOCIATED(t0_state%wf))
644  cpassert(ASSOCIATED(t1_state%wf))
645  nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
646  CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
647 
648  my_orthogonal_wf = .true.
649  t0 = 0.0_dp
650  t1 = t1_state%dt
651  t2 = t1 + dt
652  CALL qs_rho_get(rho, rho_ao=rho_ao)
653  DO ispin = 1, SIZE(mos)
654  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, &
655  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  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  beta=(t1 - t2)/(t1 - t0), matrix_b=t0_state%wf(ispin))
664  CALL reorthogonalize_vectors(qs_env, &
665  v_matrix=mo_coeff, &
666  n_col=nmo)
667  CALL calculate_density_matrix(mo_set=mos(ispin), &
668  density_matrix=rho_ao(ispin)%matrix)
669  END DO
670  CALL qs_rho_update_rho(rho, qs_env=qs_env)
671 
672  CALL qs_ks_did_change(qs_env%ks_env, &
673  rho_changed=.true.)
675  t0_state => wfi_get_snapshot(wf_history, wf_index=2)
676  t1_state => wfi_get_snapshot(wf_history, wf_index=1)
677  cpassert(ASSOCIATED(t0_state))
678  cpassert(ASSOCIATED(t1_state))
679  IF (do_kpoints) THEN
680  cpassert(ASSOCIATED(t0_state%rho_ao_kp))
681  cpassert(ASSOCIATED(t1_state%rho_ao_kp))
682  ELSE
683  cpassert(ASSOCIATED(t0_state%rho_ao))
684  cpassert(ASSOCIATED(t1_state%rho_ao))
685  END IF
686  nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
687  CALL wfi_set_history_variables(qs_env=qs_env, nvec=nvec)
688 
689  t0 = 0.0_dp
690  t1 = t1_state%dt
691  t2 = t1 + dt
692  IF (do_kpoints) THEN
693  CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
694  DO ispin = 1, SIZE(rho_ao_kp, 1)
695  DO img = 1, SIZE(rho_ao_kp, 2)
696  IF (img > SIZE(t0_state%rho_ao_kp, 2) .OR. &
697  img > SIZE(t1_state%rho_ao_kp, 2)) THEN
698  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  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  alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0))
704  END IF
705  END DO
706  END DO
707  ELSE
708  CALL qs_rho_get(rho, rho_ao=rho_ao)
709  DO ispin = 1, SIZE(rho_ao)
710  CALL dbcsr_add(rho_ao(ispin)%matrix, t1_state%rho_ao(ispin)%matrix, &
711  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  alpha_scalar=1.0_dp, beta_scalar=(t1 - t2)/(t1 - t0))
714  END DO
715  END IF
716  CALL qs_rho_update_rho(rho, qs_env=qs_env)
717  CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
719  ! wf not calculated, extract with PSC renormalized?
720  ! use wf_linear?
721  cpassert(.NOT. do_kpoints)
722  t0_state => wfi_get_snapshot(wf_history, wf_index=2)
723  t1_state => wfi_get_snapshot(wf_history, wf_index=1)
724  cpassert(ASSOCIATED(t0_state))
725  cpassert(ASSOCIATED(t1_state))
726  cpassert(ASSOCIATED(t0_state%wf))
727  cpassert(ASSOCIATED(t1_state%wf))
728  IF (wf_history%store_overlap) THEN
729  cpassert(ASSOCIATED(t0_state%overlap))
730  cpassert(ASSOCIATED(t1_state%overlap))
731  END IF
732  nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
733  IF (nvec >= wf_history%memory_depth) THEN
734  IF ((qs_env%scf_control%max_scf_hist .NE. 0) .AND. (qs_env%scf_control%eps_scf_hist .NE. 0)) THEN
735  qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
736  qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
737  qs_env%scf_control%outer_scf%have_scf = .false.
738  ELSE IF (qs_env%scf_control%max_scf_hist .NE. 0) THEN
739  qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
740  qs_env%scf_control%outer_scf%have_scf = .false.
741  ELSE IF (qs_env%scf_control%eps_scf_hist .NE. 0) THEN
742  qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
743  END IF
744  END IF
745 
746  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  CALL qs_rho_get(rho, rho_ao=rho_ao)
750  DO ispin = 1, SIZE(mos)
751  NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
752  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  matrix_struct=matrix_struct)
755  CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, &
756  nrow_global=k, ncol_global=k)
757  CALL cp_fm_create(csc, matrix_struct_new)
758  CALL cp_fm_struct_release(matrix_struct_new)
759 
760  IF (use_overlap) THEN
761  CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), mo_coeff, k)
762  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  t1_state%wf(ispin), 0.0_dp, csc)
766  END IF
767  CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, mo_coeff)
768  CALL cp_fm_release(csc)
769  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  n_col=k)
773  CALL calculate_density_matrix(mo_set=mos(ispin), &
774  density_matrix=rho_ao(ispin)%matrix)
775  END DO
776  CALL qs_rho_update_rho(rho, qs_env=qs_env)
777  CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
778 
779  CASE (wfi_ps_method_nr)
780  cpassert(.NOT. do_kpoints)
781  ! figure out the actual number of vectors to use in the extrapolation:
782  nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
783  cpassert(nvec .GT. 0)
784  IF (nvec >= wf_history%memory_depth) THEN
785  IF ((qs_env%scf_control%max_scf_hist .NE. 0) .AND. (qs_env%scf_control%eps_scf_hist .NE. 0)) THEN
786  qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
787  qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
788  qs_env%scf_control%outer_scf%have_scf = .false.
789  ELSE IF (qs_env%scf_control%max_scf_hist .NE. 0) THEN
790  qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
791  qs_env%scf_control%outer_scf%have_scf = .false.
792  ELSE IF (qs_env%scf_control%eps_scf_hist .NE. 0) THEN
793  qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
794  END IF
795  END IF
796 
797  my_orthogonal_wf = .true.
798  DO ispin = 1, SIZE(mos)
799  NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
800  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  matrix_struct=matrix_struct)
803  CALL cp_fm_create(fm_tmp, matrix_struct)
804  CALL cp_fm_struct_create(matrix_struct_new, template_fmstruct=matrix_struct, &
805  nrow_global=k, ncol_global=k)
806  CALL cp_fm_create(csc, matrix_struct_new)
807  CALL cp_fm_struct_release(matrix_struct_new)
808  ! first the most recent
809  t1_state => wfi_get_snapshot(wf_history, wf_index=1)
810  CALL cp_fm_to_fm(t1_state%wf(ispin), mo_coeff)
811  alpha = nvec
812  CALL cp_fm_scale(alpha, mo_coeff)
813  CALL qs_rho_get(rho, rho_ao=rho_ao)
814  DO i = 2, nvec
815  t0_state => wfi_get_snapshot(wf_history, wf_index=i)
816  IF (use_overlap) THEN
817  CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
818  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  t1_state%wf(ispin), 0.0_dp, csc)
822  END IF
823  CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, t0_state%wf(ispin), csc, 0.0_dp, fm_tmp)
824  alpha = -1.0_dp*alpha*real(nvec - i + 1, dp)/real(i, dp)
825  CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp)
826  END DO
827 
828  CALL cp_fm_release(csc)
829  CALL cp_fm_release(fm_tmp)
830  CALL reorthogonalize_vectors(qs_env, &
831  v_matrix=mo_coeff, &
832  n_col=k)
833  CALL calculate_density_matrix(mo_set=mos(ispin), &
834  density_matrix=rho_ao(ispin)%matrix)
835  END DO
836  CALL qs_rho_update_rho(rho, qs_env=qs_env)
837  CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
838 
839  CASE (wfi_aspc_nr)
840  cpassert(.NOT. do_kpoints)
841  CALL cite_reference(kolafa2004)
842  CALL cite_reference(kuhne2007)
843  ! figure out the actual number of vectors to use in the extrapolation:
844  nvec = min(wf_history%memory_depth, wf_history%snapshot_count)
845  cpassert(nvec .GT. 0)
846  IF (nvec >= wf_history%memory_depth) THEN
847  IF ((qs_env%scf_control%max_scf_hist .NE. 0) .AND. &
848  (qs_env%scf_control%eps_scf_hist .NE. 0)) THEN
849  qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
850  qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
851  qs_env%scf_control%outer_scf%have_scf = .false.
852  ELSE IF (qs_env%scf_control%max_scf_hist .NE. 0) THEN
853  qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
854  qs_env%scf_control%outer_scf%have_scf = .false.
855  ELSE IF (qs_env%scf_control%eps_scf_hist .NE. 0) THEN
856  qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
857  END IF
858  END IF
859 
860  my_orthogonal_wf = .true.
861  CALL qs_rho_get(rho, rho_ao=rho_ao)
862  DO ispin = 1, SIZE(mos)
863  NULLIFY (mo_coeff, matrix_struct, matrix_struct_new)
864  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  matrix_struct=matrix_struct)
869  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  ncol_global=k)
874  CALL cp_fm_create(csc, matrix_struct_new)
875  CALL cp_fm_struct_release(matrix_struct_new)
876  ! first the most recent
877  t1_state => wfi_get_snapshot(wf_history, &
878  wf_index=1)
879  CALL cp_fm_to_fm(t1_state%wf(ispin), mo_coeff)
880  alpha = real(4*nvec - 2, kind=dp)/real(nvec + 1, kind=dp)
881  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  "Parameters for the always stable predictor-corrector (ASPC) method:", &
884  "ASPC order: ", max(nvec - 2, 0), &
885  "B(", 1, ") = ", alpha
886  END IF
887  CALL cp_fm_scale(alpha, mo_coeff)
888 
889  DO i = 2, nvec
890  t0_state => wfi_get_snapshot(wf_history, wf_index=i)
891  IF (use_overlap) THEN
892  CALL cp_dbcsr_sm_fm_multiply(t0_state%overlap, t1_state%wf(ispin), fm_tmp, k)
893  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  t1_state%wf(ispin), 0.0_dp, csc)
897  END IF
898  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  binomial(2*nvec, nvec - i)/binomial(2*nvec - 2, nvec - 1)
901  IF ((io_unit > 0) .AND. (print_level > low_print_level)) THEN
902  WRITE (unit=io_unit, fmt="(T3,A2,I0,A4,F10.6)") &
903  "B(", i, ") = ", alpha
904  END IF
905  CALL cp_fm_scale_and_add(1.0_dp, mo_coeff, alpha, fm_tmp)
906  END DO
907  CALL cp_fm_release(csc)
908  CALL cp_fm_release(fm_tmp)
909  CALL reorthogonalize_vectors(qs_env, &
910  v_matrix=mo_coeff, &
911  n_col=k)
912  CALL calculate_density_matrix(mo_set=mos(ispin), &
913  density_matrix=rho_ao(ispin)%matrix)
914  END DO
915  CALL qs_rho_update_rho(rho, qs_env=qs_env)
916  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  trim(adjustl(cp_to_string(wf_history%interpolation_method_nr))))
922  END SELECT
923  IF (PRESENT(orthogonal_wf)) orthogonal_wf = my_orthogonal_wf
924  CALL cp_print_key_finished_output(io_unit, logger, qs_env%input, &
925  "DFT%SCF%PRINT%PROGRAM_RUN_INFO")
926  CALL timestop(handle)
927  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  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  IF (nvec >= qs_env%wf_history%memory_depth) THEN
943  IF ((qs_env%scf_control%max_scf_hist .NE. 0) .AND. (qs_env%scf_control%eps_scf_hist .NE. 0)) THEN
944  qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
945  qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
946  qs_env%scf_control%outer_scf%have_scf = .false.
947  ELSE IF (qs_env%scf_control%max_scf_hist .NE. 0) THEN
948  qs_env%scf_control%max_scf = qs_env%scf_control%max_scf_hist
949  qs_env%scf_control%outer_scf%have_scf = .false.
950  ELSE IF (qs_env%scf_control%eps_scf_hist .NE. 0) THEN
951  qs_env%scf_control%eps_scf = qs_env%scf_control%eps_scf_hist
952  qs_env%scf_control%outer_scf%eps_scf = qs_env%scf_control%eps_scf_hist
953  END IF
954  END IF
955 
956  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  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  cpassert(ASSOCIATED(wf_history))
973  cpassert(wf_history%ref_count > 0)
974  cpassert(ASSOCIATED(qs_env))
975 
976  wf_history%snapshot_count = wf_history%snapshot_count + 1
977  IF (wf_history%memory_depth > 0) THEN
978  wf_history%last_state_index = modulo(wf_history%snapshot_count, &
979  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  qs_env=qs_env, dt=dt)
983  END IF
984  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  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  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  NULLIFY (scf_env, scf_control, maxao_maxmo_fm_pool, matrix_s, mpools, dft_control)
1014  CALL timeset(routinen, handle)
1015 
1016  cpassert(ASSOCIATED(qs_env))
1017 
1018  CALL cp_fm_get_info(v_matrix, ncol_global=my_n_col)
1019  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  dft_control=dft_control)
1025  CALL mpools_get(mpools, maxao_maxmo_fm_pool=maxao_maxmo_fm_pool)
1026  IF (ASSOCIATED(scf_env)) THEN
1027  ortho_contains_cholesky = (scf_env%method /= ot_method_nr) .AND. &
1028  (scf_env%cholesky_method > 0) .AND. &
1029  ASSOCIATED(scf_env%ortho)
1030  ELSE
1031  ortho_contains_cholesky = .false.
1032  END IF
1033 
1034  CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
1035  smearing_is_used = .false.
1036  IF (dft_control%smear) THEN
1037  smearing_is_used = .true.
1038  END IF
1039 
1040  IF (has_unit_metric) THEN
1041  CALL make_basis_simple(v_matrix, my_n_col)
1042  ELSE IF (smearing_is_used) THEN
1043  CALL make_basis_lowdin(vmatrix=v_matrix, ncol=my_n_col, &
1044  matrix_s=matrix_s(1)%matrix)
1045  ELSE IF (ortho_contains_cholesky) THEN
1046  CALL make_basis_cholesky(vmatrix=v_matrix, ncol=my_n_col, &
1047  ortho=scf_env%ortho)
1048  ELSE
1049  CALL make_basis_sm(v_matrix, my_n_col, matrix_s(1)%matrix)
1050  END IF
1051  CALL timestop(handle)
1052  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  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  NULLIFY (dft_control, wf_history)
1073 
1074  CALL timeset(routinen, handle)
1075  logger => cp_get_default_logger()
1076  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  extension=".scfLog")
1079 
1080  cpassert(ASSOCIATED(qs_env))
1081  cpassert(ASSOCIATED(qs_env%wf_history))
1082  cpassert(qs_env%wf_history%ref_count > 0)
1083  CALL get_qs_env(qs_env, dft_control=dft_control)
1084 
1085  SELECT CASE (qs_env%wf_history%interpolation_method_nr)
1089  ! do nothing
1092  wfi_aspc_nr)
1093  IF (qs_env%wf_history%snapshot_count .GE. 2) THEN
1094  IF (debug_this_module .AND. io_unit > 0) &
1095  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  has_unit_metric=qs_env%has_unit_metric)
1100  CALL set_qs_env(qs_env=qs_env, &
1101  wf_history=wf_history)
1102  CALL wfi_release(wf_history)
1103  CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
1104  END IF
1105  CASE DEFAULT
1106  cpabort("Unknown extrapolation method.")
1107  END SELECT
1108  CALL timestop(handle)
1109 
1110  END SUBROUTINE wfi_purge_history
1111 
1112 END MODULE qs_wf_history_methods
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
integer, save, public vandevondele2005a
Definition: bibliography.F:43
integer, save, public kuhne2007
Definition: bibliography.F:43
integer, save, public kolafa2004
Definition: bibliography.F:43
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
basic linear algebra operations for full matrices
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
pool for for elements that are retained and released
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
integer, parameter, public low_print_level
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public wfi_frozen_method_nr
integer, parameter, public wfi_linear_wf_method_nr
integer, parameter, public wfi_linear_p_method_nr
integer, parameter, public wfi_linear_ps_method_nr
integer, parameter, public wfi_use_prev_rho_r_method_nr
integer, parameter, public wfi_use_guess_method_nr
integer, parameter, public wfi_use_prev_wf_method_nr
integer, parameter, public wfi_ps_method_nr
integer, parameter, public wfi_aspc_nr
integer, parameter, public wfi_use_prev_p_method_nr
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
elemental real(kind=dp) function, public binomial(n, k)
The binomial coefficient n over k for 0 <= k <= n is calculated, otherwise zero is returned.
Definition: mathlib.F:205
basic linear algebra operations for full matrixes
container for various plainwaves related things
Definition: pw_env_types.F:14
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
collects routines that calculate density matrices
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Set the QUICKSTEP environment.
subroutine, public qs_ks_did_change(ks_env, s_mstruct_changed, rho_changed, potential_changed, full_reset)
tells that some of the things relevant to the ks calculation did change. has to be called when change...
Definition: qs_ks_types.F:919
wrapper for the pools of matrixes
subroutine, public mpools_get(mpools, ao_mo_fm_pools, ao_ao_fm_pools, mo_mo_fm_pools, ao_mosub_fm_pools, mosub_mosub_fm_pools, maxao_maxmo_fm_pool, maxao_maxao_fm_pool, maxmo_maxmo_fm_pool)
returns various attributes of the mpools (notably the pools contained in it)
collects routines that perform operations directly related to MOs
Definition: qs_mo_methods.F:14
subroutine, public make_basis_simple(vmatrix, ncol)
given a set of vectors, return an orthogonal (C^T C == 1) set spanning the same space (notice,...
subroutine, public make_basis_lowdin(vmatrix, ncol, matrix_s)
return a set of S orthonormal vectors (C^T S C == 1) where a Loedwin transformation is applied to kee...
subroutine, public make_basis_cholesky(vmatrix, ncol, ortho)
return a set of S orthonormal vectors (C^T S C == 1) where the cholesky decomposed form of S is passe...
subroutine, public make_basis_sm(vmatrix, ncol, matrix_s)
returns an S-orthonormal basis v (v^T S v ==1)
Definition: qs_mo_methods.F:87
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
methods of the rho structure (defined in qs_rho_types)
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
subroutine, public qs_rho_set(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
...
Definition: qs_rho_types.F:308
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
module that contains the definitions of the scf types
Definition: qs_scf_types.F:14
integer, parameter, public ot_method_nr
Definition: qs_scf_types.F:51
Storage of past states of the qs_env. Methods to interpolate (or actually normally extrapolate) the n...
subroutine, public wfi_create_for_kp(wf_history)
...
subroutine, public wfi_purge_history(qs_env)
purges wf_history retaining only the latest snapshot
character(len=30) function, public wfi_get_method_label(method_nr)
returns a string describing the interpolation method
subroutine, public reorthogonalize_vectors(qs_env, v_matrix, n_col)
reorthogonalizes the mos
subroutine, public wfi_update(wf_history, qs_env, dt)
updates the snapshot buffer, taking a new snapshot
subroutine, public wfi_extrapolate(wf_history, qs_env, dt, extrapolation_method_nr, orthogonal_wf)
calculates the new starting state for the scf for the next wf optimization
subroutine, public wfi_create(wf_history, interpolation_method_nr, extrapolation_order, has_unit_metric)
...
interpolate the wavefunctions to speed up the convergence when doing MD
type(qs_wf_snapshot_type) function, pointer, public wfi_get_snapshot(wf_history, wf_index)
returns a snapshot, the first being the latest snapshot
subroutine, public wfi_release(wf_history)
releases a wf_history of a wavefunction (see doc/ReferenceCounting.html)
parameters that control an scf iteration