1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Types and set_get for real time propagation
8!>        depending on runtype and diagonalization method different
9!>        matrices are allocated
10!>        exp_H_old, exp_H_new, mos_new, mos_old contain always
11!>        real and imaginary parts of the matrices
12!>        odd index = real part (alpha, beta spin)
13!>        even index= imaginary part (alpha, beta spin)
14!> \par History
15!>      02.2014 switched to dbcsr matrices [Samuel Andermatt]
16!> \author Florian Schiffmann 02.09
17! **************************************************************************************************
18
19MODULE rt_propagation_types
20
21   USE bibliography,                    ONLY: Kunert2003,&
22                                              cite_reference
23   USE cp_control_types,                ONLY: dft_control_type,&
24                                              rtp_control_type
25   USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set,&
26                                              dbcsr_deallocate_matrix_set
27   USE cp_fm_pool_types,                ONLY: cp_fm_pool_p_type,&
28                                              fm_pool_get_el_struct
29   USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
30                                              cp_fm_struct_get,&
31                                              cp_fm_struct_release,&
32                                              cp_fm_struct_type
33   USE cp_fm_types,                     ONLY: cp_fm_create,&
34                                              cp_fm_p_type,&
35                                              cp_fm_release
36   USE cp_fm_vect,                      ONLY: cp_fm_vect_dealloc
37   USE cp_log_handling,                 ONLY: cp_to_string
38   USE dbcsr_api,                       ONLY: dbcsr_create,&
39                                              dbcsr_deallocate_matrix,&
40                                              dbcsr_init_p,&
41                                              dbcsr_p_type,&
42                                              dbcsr_type
43   USE kinds,                           ONLY: dp
44   USE qs_matrix_pools,                 ONLY: mpools_get,&
45                                              qs_matrix_pools_type
46   USE qs_mo_types,                     ONLY: get_mo_set,&
47                                              mo_set_p_type
48#include "./base/base_uses.f90"
49
50   IMPLICIT NONE
51
52   PRIVATE
53
54   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_types'
55
56   TYPE rtp_rho_type
57      TYPE(dbcsr_p_type), POINTER, DIMENSION(:) :: new
58      TYPE(dbcsr_p_type), POINTER, DIMENSION(:) :: old
59      TYPE(dbcsr_p_type), POINTER, DIMENSION(:) :: next
60   END TYPE rtp_rho_type
61
62   TYPE rtp_history_type
63      TYPE(dbcsr_p_type), POINTER, DIMENSION(:, :)               :: rho_history
64      TYPE(dbcsr_p_type), POINTER, DIMENSION(:)                 :: s_history
65      TYPE(cp_fm_p_type), POINTER, DIMENSION(:, :)                :: mo_history
66   END TYPE rtp_history_type
67
68   TYPE rtp_mos_type
69      TYPE(cp_fm_p_type), POINTER, DIMENSION(:)     :: new
70      TYPE(cp_fm_p_type), POINTER, DIMENSION(:)     :: old
71      TYPE(cp_fm_p_type), POINTER, DIMENSION(:)     :: next
72      TYPE(cp_fm_p_type), POINTER, DIMENSION(:)     :: admm
73   END TYPE rtp_mos_type
74
75   TYPE rt_prop_type
76      TYPE(dbcsr_p_type), POINTER, DIMENSION(:)  :: exp_H_old
77      TYPE(dbcsr_p_type), POINTER, DIMENSION(:)  :: exp_H_new
78      TYPE(dbcsr_p_type), POINTER, DIMENSION(:)  :: H_last_iter
79      TYPE(dbcsr_p_type), POINTER, DIMENSION(:)  :: propagator_matrix
80      TYPE(dbcsr_type), POINTER                 :: S_inv
81      TYPE(dbcsr_type), POINTER                 :: S_half
82      TYPE(dbcsr_type), POINTER                 :: S_minus_half
83      TYPE(dbcsr_type), POINTER                 :: B_mat
84      TYPE(dbcsr_p_type), POINTER, DIMENSION(:)  :: C_mat
85      TYPE(dbcsr_p_type), POINTER, DIMENSION(:)  :: S_der
86      TYPE(dbcsr_p_type), POINTER, DIMENSION(:)  :: SinvH
87      TYPE(dbcsr_p_type), POINTER, DIMENSION(:) :: SinvB
88      TYPE(rtp_rho_type), POINTER                  :: rho
89      TYPE(rtp_mos_type), POINTER                  :: mos
90      REAL(KIND=dp)                               :: energy_old
91      REAL(KIND=dp)                               :: energy_new
92      REAL(KIND=dp)                               :: dt
93      REAL(KIND=dp)                               :: delta_iter
94      REAL(KIND=dp)                               :: delta_iter_old
95      REAL(KIND=dp)                               :: filter_eps
96      REAL(KIND=dp)                               :: filter_eps_small
97      REAL(KIND=dp)                               :: mixing_factor
98      LOGICAL                                     :: mixing
99      LOGICAL                                     :: do_hfx
100      LOGICAL                                     :: magnetic
101      INTEGER, DIMENSION(:, :), ALLOCATABLE          :: orders
102      INTEGER                                     :: nsteps, istep, i_start, max_steps
103      INTEGER                                     :: iter
104      INTEGER                                     :: narn_old
105      LOGICAL                                     :: converged
106      LOGICAL                                     :: matrix_update
107      LOGICAL                                     :: write_restart
108      TYPE(rtp_history_type), POINTER              :: history
109      TYPE(cp_fm_struct_type), POINTER            :: ao_ao_fmstruct
110      INTEGER                                     :: lanzcos_max_iter
111      REAL(KIND=dp)                               :: lanzcos_threshold
112      INTEGER                                     :: newton_schulz_order
113      LOGICAL                                     :: linear_scaling
114   END TYPE rt_prop_type
115
116! *** Public data types ***
117
118   PUBLIC :: rt_prop_type
119
120! *** Public subroutines ***
121
122   PUBLIC :: rt_prop_create, &
123             get_rtp, &
124             rt_prop_release, &
125             rtp_history_create
126CONTAINS
127
128! **************************************************************************************************
129!> \brief ...
130!> \param rtp ...
131!> \param mos ...
132!> \param mpools ...
133!> \param dft_control ...
134!> \param template ...
135!> \param linear_scaling ...
136!> \param mos_aux ...
137! **************************************************************************************************
138   SUBROUTINE rt_prop_create(rtp, mos, mpools, dft_control, template, linear_scaling, mos_aux)
139
140      TYPE(rt_prop_type), POINTER                        :: rtp
141      TYPE(mo_set_p_type), DIMENSION(:), POINTER         :: mos
142      TYPE(qs_matrix_pools_type), POINTER                :: mpools
143      TYPE(dft_control_type), POINTER                    :: dft_control
144      TYPE(dbcsr_type), POINTER                          :: template
145      LOGICAL                                            :: linear_scaling
146      TYPE(mo_set_p_type), DIMENSION(:), POINTER         :: mos_aux
147
148      CHARACTER(len=*), PARAMETER :: routineN = 'rt_prop_create', routineP = moduleN//':'//routineN
149
150      INTEGER                                            :: i, j, nao, nrow_block, nspin
151      TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER     :: ao_mo_fm_pools
152      TYPE(cp_fm_struct_type), POINTER                   :: ao_mo_fmstruct
153      TYPE(rtp_control_type), POINTER                    :: rtp_control
154
155      CALL cite_reference(Kunert2003)
156
157      NULLIFY (rtp_control)
158
159      rtp_control => dft_control%rtp_control
160
161      nspin = dft_control%nspins
162
163      NULLIFY (rtp%mos, rtp%rho)
164      rtp%linear_scaling = linear_scaling
165
166      IF (rtp%linear_scaling) THEN
167         ALLOCATE (rtp%rho)
168         NULLIFY (rtp%rho%old)
169         CALL dbcsr_allocate_matrix_set(rtp%rho%old, 2*nspin)
170         NULLIFY (rtp%rho%next)
171         CALL dbcsr_allocate_matrix_set(rtp%rho%next, 2*nspin)
172         NULLIFY (rtp%rho%new)
173         CALL dbcsr_allocate_matrix_set(rtp%rho%new, 2*nspin)
174         DO i = 1, 2*nspin
175            CALL dbcsr_init_p(rtp%rho%old(i)%matrix)
176            CALL dbcsr_create(rtp%rho%old(i)%matrix, template=template, matrix_type="N")
177            CALL dbcsr_init_p(rtp%rho%next(i)%matrix)
178            CALL dbcsr_create(rtp%rho%next(i)%matrix, template=template, matrix_type="N")
179            CALL dbcsr_init_p(rtp%rho%new(i)%matrix)
180            CALL dbcsr_create(rtp%rho%new(i)%matrix, template=template, matrix_type="N")
181         END DO
182      ELSE
183         CALL mpools_get(mpools, ao_mo_fm_pools=ao_mo_fm_pools)
184
185         ao_mo_fmstruct => fm_pool_get_el_struct(ao_mo_fm_pools(1)%pool)
186         CALL cp_fm_struct_get(ao_mo_fmstruct, nrow_block=nrow_block)
187         CALL get_mo_set(mos(1)%mo_set, nao=nao)
188
189         CALL cp_fm_struct_create(fmstruct=rtp%ao_ao_fmstruct, &
190                                  nrow_block=nrow_block, ncol_block=nrow_block, &
191                                  nrow_global=nao, ncol_global=nao, &
192                                  template_fmstruct=ao_mo_fmstruct)
193         ALLOCATE (rtp%mos)
194         ALLOCATE (rtp%mos%old(2*nspin))
195         ALLOCATE (rtp%mos%new(2*nspin))
196         ALLOCATE (rtp%mos%next(2*nspin))
197         NULLIFY (rtp%mos%admm)
198         IF (dft_control%do_admm) THEN
199            ALLOCATE (rtp%mos%admm(2*nspin))
200         END IF
201         DO i = 1, nspin
202            DO j = 1, 2
203               NULLIFY (rtp%mos%old(2*(i - 1) + j)%matrix)
204               NULLIFY (rtp%mos%new(2*(i - 1) + j)%matrix)
205               NULLIFY (rtp%mos%next(2*(i - 1) + j)%matrix)
206               CALL cp_fm_create(rtp%mos%old(2*(i - 1) + j)%matrix, &
207                                 matrix_struct=mos(i)%mo_set%mo_coeff%matrix_struct, &
208                                 name="mos_old"//TRIM(ADJUSTL(cp_to_string(2*(i - 1) + j))))
209               CALL cp_fm_create(rtp%mos%new(2*(i - 1) + j)%matrix, &
210                                 matrix_struct=mos(i)%mo_set%mo_coeff%matrix_struct, &
211                                 name="mos_new"//TRIM(ADJUSTL(cp_to_string(2*(i - 1) + j))))
212               CALL cp_fm_create(rtp%mos%next(2*(i - 1) + j)%matrix, &
213                                 matrix_struct=mos(i)%mo_set%mo_coeff%matrix_struct, &
214                                 name="mos_next"//TRIM(ADJUSTL(cp_to_string(2*(i - 1) + j))))
215               IF (dft_control%do_admm) THEN
216                  NULLIFY (rtp%mos%admm(2*(i - 1) + j)%matrix)
217                  CALL cp_fm_create(rtp%mos%admm(2*(i - 1) + j)%matrix, &
218                                    matrix_struct=mos_aux(i)%mo_set%mo_coeff%matrix_struct, &
219                                    name="mos_admm"//TRIM(ADJUSTL(cp_to_string(2*(i - 1) + j))))
220               END IF
221            END DO
222         END DO
223      END IF
224
225      NULLIFY (rtp%exp_H_old)
226      NULLIFY (rtp%exp_H_new)
227      NULLIFY (rtp%H_last_iter)
228      NULLIFY (rtp%propagator_matrix)
229      CALL dbcsr_allocate_matrix_set(rtp%exp_H_old, 2*nspin)
230      CALL dbcsr_allocate_matrix_set(rtp%exp_H_new, 2*nspin)
231      CALL dbcsr_allocate_matrix_set(rtp%H_last_iter, 2*nspin)
232      CALL dbcsr_allocate_matrix_set(rtp%propagator_matrix, 2*nspin)
233      DO i = 1, 2*nspin
234         CALL dbcsr_init_p(rtp%exp_H_old(i)%matrix)
235         CALL dbcsr_create(rtp%exp_H_old(i)%matrix, template=template, matrix_type="N")
236         CALL dbcsr_init_p(rtp%exp_H_new(i)%matrix)
237         CALL dbcsr_create(rtp%exp_H_new(i)%matrix, template=template, matrix_type="N")
238         CALL dbcsr_init_p(rtp%H_last_iter(i)%matrix)
239         CALL dbcsr_create(rtp%H_last_iter(i)%matrix, template=template, matrix_type="N")
240         CALL dbcsr_init_p(rtp%propagator_matrix(i)%matrix)
241         CALL dbcsr_create(rtp%propagator_matrix(i)%matrix, template=template, matrix_type="N")
242      END DO
243      NULLIFY (rtp%S_inv)
244      ALLOCATE (rtp%S_inv)
245      CALL dbcsr_create(rtp%S_inv, template=template, matrix_type="S")
246      NULLIFY (rtp%S_half)
247      ALLOCATE (rtp%S_half)
248      CALL dbcsr_create(rtp%S_half, template=template, matrix_type="S")
249      NULLIFY (rtp%S_minus_half)
250      ALLOCATE (rtp%S_minus_half)
251      CALL dbcsr_create(rtp%S_minus_half, template=template, matrix_type="S")
252      NULLIFY (rtp%B_mat)
253      NULLIFY (rtp%C_mat)
254      NULLIFY (rtp%S_der)
255      NULLIFY (rtp%SinvH)
256      NULLIFY (rtp%SinvB)
257      IF (.NOT. rtp_control%fixed_ions) THEN
258         ALLOCATE (rtp%B_mat)
259         CALL dbcsr_create(rtp%B_mat, template=template, matrix_type="N")
260         CALL dbcsr_allocate_matrix_set(rtp%C_mat, 3)
261         CALL dbcsr_allocate_matrix_set(rtp%S_der, 9)
262         CALL dbcsr_allocate_matrix_set(rtp%SinvH, nspin)
263         CALL dbcsr_allocate_matrix_set(rtp%SinvB, nspin)
264         DO i = 1, nspin
265            CALL dbcsr_init_p(rtp%SinvH(i)%matrix)
266            CALL dbcsr_create(rtp%SinvH(i)%matrix, template=template, matrix_type="N")
267            CALL dbcsr_init_p(rtp%SinvB(i)%matrix)
268            CALL dbcsr_create(rtp%SinvB(i)%matrix, template=template, matrix_type="N")
269         END DO
270         DO i = 1, 3
271            CALL dbcsr_init_p(rtp%C_mat(i)%matrix)
272            CALL dbcsr_create(rtp%C_mat(i)%matrix, template=template, matrix_type="N")
273         END DO
274         DO i = 1, 9
275            CALL dbcsr_init_p(rtp%S_der(i)%matrix)
276            CALL dbcsr_create(rtp%S_der(i)%matrix, template=template, matrix_type="N")
277         END DO
278      END IF
279      ALLOCATE (rtp%orders(2, nspin))
280      rtp_control%converged = .FALSE.
281      rtp%matrix_update = .TRUE.
282      rtp%narn_old = 0
283      rtp%istep = 0
284      rtp%iter = 0
285      rtp%do_hfx = .FALSE.
286      rtp%magnetic = .FALSE.
287
288   END SUBROUTINE rt_prop_create
289
290! **************************************************************************************************
291!> \brief ...
292!> \param rtp ...
293!> \param exp_H_old ...
294!> \param exp_H_new ...
295!> \param H_last_iter ...
296!> \param rho_old ...
297!> \param rho_next ...
298!> \param rho_new ...
299!> \param mos ...
300!> \param mos_new ...
301!> \param mos_old ...
302!> \param mos_next ...
303!> \param S_inv ...
304!> \param S_half ...
305!> \param S_minus_half ...
306!> \param B_mat ...
307!> \param C_mat ...
308!> \param propagator_matrix ...
309!> \param mixing ...
310!> \param mixing_factor ...
311!> \param S_der ...
312!> \param dt ...
313!> \param nsteps ...
314!> \param SinvH ...
315!> \param SinvB ...
316!> \param admm_mos ...
317! **************************************************************************************************
318   SUBROUTINE get_rtp(rtp, exp_H_old, exp_H_new, H_last_iter, rho_old, rho_next, rho_new, mos, mos_new, mos_old, mos_next, &
319                      S_inv, S_half, S_minus_half, B_mat, C_mat, propagator_matrix, mixing, mixing_factor, &
320                      S_der, dt, nsteps, SinvH, SinvB, admm_mos)
321
322      TYPE(rt_prop_type), POINTER                        :: rtp
323      TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
324         POINTER                                         :: exp_H_old, exp_H_new, H_last_iter, &
325                                                            rho_old, rho_next, rho_new
326      TYPE(rtp_mos_type), OPTIONAL, POINTER              :: mos
327      TYPE(cp_fm_p_type), DIMENSION(:), OPTIONAL, &
328         POINTER                                         :: mos_new, mos_old, mos_next
329      TYPE(dbcsr_type), OPTIONAL, POINTER                :: S_inv, S_half, S_minus_half, B_mat
330      TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
331         POINTER                                         :: C_mat, propagator_matrix
332      LOGICAL, OPTIONAL                                  :: mixing
333      REAL(dp), INTENT(out), OPTIONAL                    :: mixing_factor
334      TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
335         POINTER                                         :: S_der
336      REAL(dp), INTENT(out), OPTIONAL                    :: dt
337      INTEGER, INTENT(out), OPTIONAL                     :: nsteps
338      TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
339         POINTER                                         :: SinvH, SinvB
340      TYPE(cp_fm_p_type), DIMENSION(:), OPTIONAL, &
341         POINTER                                         :: admm_mos
342
343      CHARACTER(len=*), PARAMETER :: routineN = 'get_rtp', routineP = moduleN//':'//routineN
344
345      CPASSERT(ASSOCIATED(rtp))
346      IF (PRESENT(exp_H_old)) exp_H_old => rtp%exp_H_old
347      IF (PRESENT(exp_H_new)) exp_H_new => rtp%exp_H_new
348      IF (PRESENT(H_last_iter)) H_last_iter => rtp%H_last_iter
349      IF (PRESENT(propagator_matrix)) propagator_matrix => rtp%propagator_matrix
350
351      IF (PRESENT(rho_old)) rho_old => rtp%rho%old
352      IF (PRESENT(rho_next)) rho_next => rtp%rho%next
353      IF (PRESENT(rho_new)) rho_new => rtp%rho%new
354      IF (PRESENT(mos)) mos => rtp%mos
355      IF (PRESENT(mos_old)) mos_old => rtp%mos%old
356      IF (PRESENT(mos_new)) mos_new => rtp%mos%new
357      IF (PRESENT(mos_next)) mos_next => rtp%mos%next
358      IF (PRESENT(admm_mos)) admm_mos => rtp%mos%admm
359
360      IF (PRESENT(S_inv)) S_inv => rtp%S_inv
361      IF (PRESENT(S_half)) S_half => rtp%S_half
362      IF (PRESENT(S_minus_half)) S_minus_half => rtp%S_minus_half
363      IF (PRESENT(B_mat)) B_mat => rtp%B_mat
364      IF (PRESENT(C_mat)) C_mat => rtp%C_mat
365      IF (PRESENT(SinvH)) SinvH => rtp%SinvH
366      IF (PRESENT(SinvB)) SinvB => rtp%SinvB
367      IF (PRESENT(S_der)) S_der => rtp%S_der
368
369      IF (PRESENT(dt)) dt = rtp%dt
370      IF (PRESENT(mixing)) mixing = rtp%mixing
371      IF (PRESENT(mixing_factor)) mixing_factor = rtp%mixing_factor
372      IF (PRESENT(nsteps)) nsteps = rtp%nsteps
373
374   END SUBROUTINE get_rtp
375
376! **************************************************************************************************
377!> \brief ...
378!> \param rtp ...
379! **************************************************************************************************
380   SUBROUTINE rt_prop_release(rtp)
381      TYPE(rt_prop_type), INTENT(inout)                  :: rtp
382
383      CHARACTER(len=*), PARAMETER :: routineN = 'rt_prop_release', &
384         routineP = moduleN//':'//routineN
385
386      CALL dbcsr_deallocate_matrix_set(rtp%exp_H_old)
387      CALL dbcsr_deallocate_matrix_set(rtp%exp_H_new)
388      CALL dbcsr_deallocate_matrix_set(rtp%H_last_iter)
389      CALL dbcsr_deallocate_matrix_set(rtp%propagator_matrix)
390      IF (ASSOCIATED(rtp%rho)) THEN
391         IF (ASSOCIATED(rtp%rho%old)) &
392            CALL dbcsr_deallocate_matrix_set(rtp%rho%old)
393         IF (ASSOCIATED(rtp%rho%next)) &
394            CALL dbcsr_deallocate_matrix_set(rtp%rho%next)
395         IF (ASSOCIATED(rtp%rho%new)) &
396            CALL dbcsr_deallocate_matrix_set(rtp%rho%new)
397         DEALLOCATE (rtp%rho)
398      ENDIF
399      IF (ASSOCIATED(rtp%mos)) THEN
400         IF (ASSOCIATED(rtp%mos%old)) &
401            CALL cp_fm_vect_dealloc(rtp%mos%old)
402         IF (ASSOCIATED(rtp%mos%new)) &
403            CALL cp_fm_vect_dealloc(rtp%mos%new)
404         IF (ASSOCIATED(rtp%mos%next)) &
405            CALL cp_fm_vect_dealloc(rtp%mos%next)
406         IF (ASSOCIATED(rtp%mos%admm)) &
407            CALL cp_fm_vect_dealloc(rtp%mos%admm)
408         DEALLOCATE (rtp%mos)
409      END IF
410      CALL dbcsr_deallocate_matrix(rtp%S_inv)
411      CALL dbcsr_deallocate_matrix(rtp%S_half)
412      CALL dbcsr_deallocate_matrix(rtp%S_minus_half)
413      IF (ASSOCIATED(rtp%B_mat)) &
414         CALL dbcsr_deallocate_matrix(rtp%B_mat)
415      IF (ASSOCIATED(rtp%C_mat)) &
416         CALL dbcsr_deallocate_matrix_set(rtp%C_mat)
417      IF (ASSOCIATED(rtp%S_der)) &
418         CALL dbcsr_deallocate_matrix_set(rtp%S_der)
419      IF (ASSOCIATED(rtp%SinvH)) &
420         CALL dbcsr_deallocate_matrix_set(rtp%SinvH)
421      IF (ASSOCIATED(rtp%SinvB)) &
422         CALL dbcsr_deallocate_matrix_set(rtp%SinvB)
423      IF (ASSOCIATED(rtp%history)) &
424         CALL rtp_history_release(rtp)
425      DEALLOCATE (rtp%orders)
426      IF (.NOT. rtp%linear_scaling) CALL cp_fm_struct_release(rtp%ao_ao_fmstruct)
427   END SUBROUTINE rt_prop_release
428
429! **************************************************************************************************
430!> \brief ...
431!> \param rtp ...
432!> \param aspc_order ...
433! **************************************************************************************************
434   SUBROUTINE rtp_history_create(rtp, aspc_order)
435      TYPE(rt_prop_type), INTENT(inout)                  :: rtp
436      INTEGER, INTENT(in)                                :: aspc_order
437
438      CHARACTER(len=*), PARAMETER :: routineN = 'rtp_history_create', &
439         routineP = moduleN//':'//routineN
440
441      INTEGER                                            :: i, j, nmat
442      TYPE(rtp_history_type), POINTER                    :: history
443
444      NULLIFY (history)
445      ALLOCATE (rtp%history)
446      history => rtp%history
447
448      NULLIFY (history%rho_history, history%mo_history, history%s_history)
449      IF (aspc_order .GT. 0) THEN
450         IF (rtp%linear_scaling) THEN
451            nmat = SIZE(rtp%rho%new)
452            CALL dbcsr_allocate_matrix_set(history%rho_history, nmat, aspc_order)
453            DO i = 1, nmat
454               DO j = 1, aspc_order
455                  CALL dbcsr_init_p(history%rho_history(i, j)%matrix)
456                  CALL dbcsr_create(history%rho_history(i, j)%matrix, &
457                                    name="rho_hist"//TRIM(ADJUSTL(cp_to_string(i))), &
458                                    template=rtp%rho%new(1)%matrix)
459               END DO
460            END DO
461         ELSE
462            nmat = SIZE(rtp%mos%old)
463            ALLOCATE (history%mo_history(nmat, aspc_order))
464            DO i = 1, nmat
465               DO j = 1, aspc_order
466                  NULLIFY (history%mo_history(i, j)%matrix)
467                  CALL cp_fm_create(history%mo_history(i, j)%matrix, &
468                                    matrix_struct=rtp%mos%new(i)%matrix%matrix_struct, &
469                                    name="mo_hist"//TRIM(ADJUSTL(cp_to_string(i))))
470               END DO
471            END DO
472            ALLOCATE (history%s_history(aspc_order))
473            DO i = 1, aspc_order
474               NULLIFY (history%s_history(i)%matrix)
475            END DO
476         END IF
477      END IF
478
479   END SUBROUTINE rtp_history_create
480
481! **************************************************************************************************
482!> \brief ...
483!> \param rtp ...
484! **************************************************************************************************
485   SUBROUTINE rtp_history_release(rtp)
486      TYPE(rt_prop_type), INTENT(inout)                  :: rtp
487
488      CHARACTER(len=*), PARAMETER :: routineN = 'rtp_history_release', &
489         routineP = moduleN//':'//routineN
490
491      INTEGER                                            :: i, j
492
493      IF (ASSOCIATED(rtp%history%rho_history)) THEN
494         CALL dbcsr_deallocate_matrix_set(rtp%history%rho_history)
495      END IF
496
497      IF (ASSOCIATED(rtp%history%mo_history)) THEN
498         DO i = 1, SIZE(rtp%history%mo_history, 1)
499            DO j = 1, SIZE(rtp%history%mo_history, 2)
500               CALL cp_fm_release(rtp%history%mo_history(i, j)%matrix)
501            END DO
502         END DO
503         DEALLOCATE (rtp%history%mo_history)
504      END IF
505      IF (ASSOCIATED(rtp%history%s_history)) THEN
506         DO i = 1, SIZE(rtp%history%s_history)
507            IF (ASSOCIATED(rtp%history%s_history(i)%matrix)) &
508               CALL dbcsr_deallocate_matrix(rtp%history%s_history(i)%matrix)
509         END DO
510         DEALLOCATE (rtp%history%s_history)
511      END IF
512      DEALLOCATE (rtp%history)
513
514   END SUBROUTINE rtp_history_release
515
516END MODULE rt_propagation_types
517