1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Routines for calculating a complex matrix exponential.
8!> \author Florian Schiffmann (02.09)
9! **************************************************************************************************
10
11MODULE rt_make_propagators
12
13   USE cp_control_types,                ONLY: rtp_control_type
14   USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
15                                              copy_fm_to_dbcsr,&
16                                              cp_dbcsr_sm_fm_multiply
17   USE cp_fm_types,                     ONLY: cp_fm_create,&
18                                              cp_fm_get_info,&
19                                              cp_fm_p_type,&
20                                              cp_fm_to_fm
21   USE cp_fm_vect,                      ONLY: cp_fm_vect_dealloc
22   USE dbcsr_api,                       ONLY: dbcsr_copy,&
23                                              dbcsr_create,&
24                                              dbcsr_deallocate_matrix,&
25                                              dbcsr_p_type,&
26                                              dbcsr_scale,&
27                                              dbcsr_type
28   USE input_constants,                 ONLY: do_etrs,&
29                                              do_pade,&
30                                              do_taylor
31   USE kinds,                           ONLY: dp
32   USE ls_matrix_exp,                   ONLY: bch_expansion_complex_propagator,&
33                                              bch_expansion_imaginary_propagator,&
34                                              cp_complex_dbcsr_gemm_3,&
35                                              taylor_full_complex_dbcsr,&
36                                              taylor_only_imaginary_dbcsr
37   USE matrix_exp,                      ONLY: arnoldi,&
38                                              exp_pade_full_complex,&
39                                              exp_pade_only_imaginary,&
40                                              taylor_full_complex,&
41                                              taylor_only_imaginary
42   USE rt_propagation_types,            ONLY: get_rtp,&
43                                              rt_prop_type
44#include "../base/base_uses.f90"
45
46   IMPLICIT NONE
47
48   PRIVATE
49
50   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_make_propagators'
51
52   PUBLIC :: propagate_exp, &
53             propagate_arnoldi, &
54             compute_exponential, &
55             compute_exponential_sparse, &
56             propagate_exp_density, &
57             propagate_bch
58
59CONTAINS
60! **************************************************************************************************
61!> \brief performs propagations if explicit matrix exponentials are used
62!>        ETRS:  exp(i*H(t+dt)*dt/2)*exp(i*H(t)*dt/2)*MOS
63!>        EM:    exp[-idt/2H(t+dt/2)*MOS
64!> \param rtp ...
65!> \param rtp_control ...
66!> \author Florian Schiffmann (02.09)
67! **************************************************************************************************
68
69   SUBROUTINE propagate_exp(rtp, rtp_control)
70
71      TYPE(rt_prop_type), POINTER                        :: rtp
72      TYPE(rtp_control_type), POINTER                    :: rtp_control
73
74      CHARACTER(len=*), PARAMETER :: routineN = 'propagate_exp', routineP = moduleN//':'//routineN
75      REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
76
77      INTEGER                                            :: handle, i, im, nmo, re
78      TYPE(cp_fm_p_type), DIMENSION(:), POINTER          :: mos_new, mos_next, mos_old
79      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: exp_H_new, exp_H_old, propagator_matrix
80
81      CALL timeset(routineN, handle)
82
83      CALL get_rtp(rtp=rtp, propagator_matrix=propagator_matrix, mos_old=mos_old, mos_new=mos_new, &
84                   mos_next=mos_next, exp_H_new=exp_H_new, exp_H_old=exp_H_old)
85
86      ! Only compute exponential if a new propagator matrix is available
87      CALL compute_exponential(exp_H_new, propagator_matrix, rtp_control, rtp)
88
89      DO i = 1, SIZE(mos_new)/2
90         re = 2*i - 1
91         im = 2*i
92
93         CALL cp_fm_get_info(mos_new(re)%matrix, ncol_global=nmo)
94         !Save some work by computing the first half of the propagation only once in case of ETRS
95         !For EM this matrix has to be the initial matrix, thus a copy is enough
96         IF (rtp%iter == 1) THEN
97            IF (rtp_control%propagator == do_etrs) THEN
98               CALL cp_dbcsr_sm_fm_multiply(exp_H_old(re)%matrix, mos_old(re)%matrix, &
99                                            mos_next(re)%matrix, nmo, alpha=one, beta=zero)
100               CALL cp_dbcsr_sm_fm_multiply(exp_H_old(im)%matrix, mos_old(im)%matrix, &
101                                            mos_next(re)%matrix, nmo, alpha=-one, beta=one)
102               CALL cp_dbcsr_sm_fm_multiply(exp_H_old(re)%matrix, mos_old(im)%matrix, &
103                                            mos_next(im)%matrix, nmo, alpha=one, beta=zero)
104               CALL cp_dbcsr_sm_fm_multiply(exp_H_old(im)%matrix, mos_old(re)%matrix, &
105                                            mos_next(im)%matrix, nmo, alpha=one, beta=one)
106            ELSE
107               CALL cp_fm_to_fm(mos_old(re)%matrix, mos_next(re)%matrix)
108               CALL cp_fm_to_fm(mos_old(im)%matrix, mos_next(im)%matrix)
109            END IF
110         END IF
111         CALL cp_dbcsr_sm_fm_multiply(exp_H_new(re)%matrix, mos_next(re)%matrix, &
112                                      mos_new(re)%matrix, nmo, alpha=one, beta=zero)
113         CALL cp_dbcsr_sm_fm_multiply(exp_H_new(im)%matrix, mos_next(im)%matrix, &
114                                      mos_new(re)%matrix, nmo, alpha=-one, beta=one)
115         CALL cp_dbcsr_sm_fm_multiply(exp_H_new(re)%matrix, mos_next(im)%matrix, &
116                                      mos_new(im)%matrix, nmo, alpha=one, beta=zero)
117         CALL cp_dbcsr_sm_fm_multiply(exp_H_new(im)%matrix, mos_next(re)%matrix, &
118                                      mos_new(im)%matrix, nmo, alpha=one, beta=one)
119      END DO
120
121      CALL timestop(handle)
122
123   END SUBROUTINE propagate_exp
124
125! **************************************************************************************************
126!> \brief Propagation of the density matrix instead of the atomic orbitals
127!>        via a matrix exponential
128!> \param rtp ...
129!> \param rtp_control ...
130!> \author Samuel Andermatt (02.2014)
131! **************************************************************************************************
132
133   SUBROUTINE propagate_exp_density(rtp, rtp_control)
134
135      TYPE(rt_prop_type), POINTER                        :: rtp
136      TYPE(rtp_control_type), POINTER                    :: rtp_control
137
138      CHARACTER(len=*), PARAMETER :: routineN = 'propagate_exp_density', &
139         routineP = moduleN//':'//routineN
140      REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
141
142      INTEGER                                            :: handle, i, im, re
143      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: exp_H_new, exp_H_old, propagator_matrix, &
144                                                            rho_new, rho_next, rho_old
145      TYPE(dbcsr_type), POINTER                          :: tmp_im, tmp_re
146
147      CALL timeset(routineN, handle)
148
149      CALL get_rtp(rtp=rtp, propagator_matrix=propagator_matrix, exp_H_new=exp_H_new, &
150                   exp_H_old=exp_H_old, rho_old=rho_old, rho_new=rho_new, rho_next=rho_next)
151
152      CALL compute_exponential_sparse(exp_H_new, propagator_matrix, rtp_control, rtp)
153
154      !I could store these matrices in the type
155      NULLIFY (tmp_re)
156      ALLOCATE (tmp_re)
157      CALL dbcsr_create(tmp_re, template=propagator_matrix(1)%matrix, matrix_type="N")
158      NULLIFY (tmp_im)
159      ALLOCATE (tmp_im)
160      CALL dbcsr_create(tmp_im, template=propagator_matrix(1)%matrix, matrix_type="N")
161
162      DO i = 1, SIZE(exp_H_new)/2
163         re = 2*i - 1
164         im = 2*i
165         !Save some work by computing the first half of the propagation only once in case of ETRS
166         !For EM this matrix has to be the initial matrix, thus a copy is enough
167         IF (rtp%iter == 1) THEN
168            IF (rtp_control%propagator == do_etrs) THEN
169               CALL cp_complex_dbcsr_gemm_3("N", "N", one, exp_H_old(re)%matrix, exp_H_old(im)%matrix, &
170                                            rho_old(re)%matrix, rho_old(im)%matrix, zero, tmp_re, tmp_im, filter_eps=rtp%filter_eps)
171               CALL cp_complex_dbcsr_gemm_3("N", "C", one, tmp_re, tmp_im, exp_H_old(re)%matrix, exp_H_old(im)%matrix, &
172                                            zero, rho_next(re)%matrix, rho_next(im)%matrix, filter_eps=rtp%filter_eps)
173            ELSE
174               CALL dbcsr_copy(rho_next(re)%matrix, rho_old(re)%matrix)
175               CALL dbcsr_copy(rho_next(im)%matrix, rho_old(im)%matrix)
176            ENDIF
177         END IF
178         CALL cp_complex_dbcsr_gemm_3("N", "N", one, exp_H_new(re)%matrix, exp_H_new(im)%matrix, &
179                                      rho_next(re)%matrix, rho_next(im)%matrix, zero, tmp_re, tmp_im, filter_eps=rtp%filter_eps)
180         CALL cp_complex_dbcsr_gemm_3("N", "C", one, tmp_re, tmp_im, exp_H_new(re)%matrix, exp_H_new(im)%matrix, &
181                                      zero, rho_new(re)%matrix, rho_new(im)%matrix, filter_eps=rtp%filter_eps)
182      END DO
183
184      CALL dbcsr_deallocate_matrix(tmp_re)
185      CALL dbcsr_deallocate_matrix(tmp_im)
186
187      CALL timestop(handle)
188
189   END SUBROUTINE propagate_exp_density
190
191! **************************************************************************************************
192!> \brief computes U_prop*MOs using arnoldi subspace algorithm
193!> \param rtp ...
194!> \param rtp_control ...
195!> \author Florian Schiffmann (02.09)
196! **************************************************************************************************
197
198   SUBROUTINE propagate_arnoldi(rtp, rtp_control)
199      TYPE(rt_prop_type), POINTER                        :: rtp
200      TYPE(rtp_control_type), POINTER                    :: rtp_control
201
202      CHARACTER(len=*), PARAMETER :: routineN = 'propagate_arnoldi', &
203         routineP = moduleN//':'//routineN
204
205      INTEGER                                            :: handle, i, im, ispin, nspin, re
206      REAL(dp)                                           :: eps_arnoldi, t
207      TYPE(cp_fm_p_type), DIMENSION(:), POINTER          :: mos_new, mos_next, mos_old, &
208                                                            propagator_matrix_fm
209      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: propagator_matrix
210
211      CALL timeset(routineN, handle)
212
213      CALL get_rtp(rtp=rtp, dt=t, mos_new=mos_new, mos_old=mos_old, &
214                   mos_next=mos_next, propagator_matrix=propagator_matrix)
215
216      nspin = SIZE(mos_new)/2
217      eps_arnoldi = rtp_control%eps_exp
218      ! except for the first step the further propagated mos_next
219      ! must be copied on mos_old so that we have the half propagated mos
220      ! ready on mos_old and only need to perform the second half propagatioon
221      IF (rtp_control%propagator == do_etrs .AND. rtp%iter == 1) THEN
222         DO i = 1, SIZE(mos_new)
223            CALL cp_fm_to_fm(mos_next(i)%matrix, mos_old(i)%matrix)
224         END DO
225      END IF
226
227      NULLIFY (propagator_matrix_fm)
228      ALLOCATE (propagator_matrix_fm(SIZE(propagator_matrix)))
229      DO i = 1, SIZE(propagator_matrix)
230         CALL cp_fm_create(propagator_matrix_fm(i)%matrix, &
231                           matrix_struct=rtp%ao_ao_fmstruct, &
232                           name="prop_fm")
233         CALL copy_dbcsr_to_fm(propagator_matrix(i)%matrix, propagator_matrix_fm(i)%matrix)
234      END DO
235
236      DO ispin = 1, nspin
237         re = ispin*2 - 1
238         im = ispin*2
239         IF (rtp_control%fixed_ions .AND. .NOT. rtp%do_hfx) THEN
240            CALL arnoldi(mos_old(re:im), mos_new(re:im), &
241                         eps_arnoldi, Him=propagator_matrix_fm(im)%matrix, &
242                         mos_next=mos_next(re:im), narn_old=rtp%narn_old)
243         ELSE
244            CALL arnoldi(mos_old(re:im), mos_new(re:im), &
245                         eps_arnoldi, Hre=propagator_matrix_fm(re)%matrix, &
246                         Him=propagator_matrix_fm(im)%matrix, mos_next=mos_next(re:im), &
247                         narn_old=rtp%narn_old)
248         END IF
249      END DO
250
251!    DO i=1,SIZE(propagator_matrix)
252!         CALL copy_fm_to_dbcsr(propagator_matrix_fm(i)%matrix,propagator_matrix(i)%matrix)
253!    END DO
254      CALL cp_fm_vect_dealloc(propagator_matrix_fm)
255
256      CALL timestop(handle)
257
258   END SUBROUTINE propagate_arnoldi
259
260! **************************************************************************************************
261!> \brief  Propagation using the Baker-Campbell-Hausdorff expansion,
262!>         currently only works for rtp
263!> \param rtp ...
264!> \param rtp_control ...
265!> \author Samuel Andermatt (02.2014)
266! **************************************************************************************************
267
268   SUBROUTINE propagate_bch(rtp, rtp_control)
269
270      TYPE(rt_prop_type), POINTER                        :: rtp
271      TYPE(rtp_control_type), POINTER                    :: rtp_control
272
273      CHARACTER(len=*), PARAMETER :: routineN = 'propagate_bch', routineP = moduleN//':'//routineN
274
275      INTEGER                                            :: handle, im, ispin, re
276      REAL(dp)                                           :: dt
277      REAL(KIND=dp)                                      :: prefac
278      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: exp_H_old, propagator_matrix, rho_new, &
279                                                            rho_next, rho_old
280
281      CALL timeset(routineN, handle)
282
283      CALL get_rtp(rtp=rtp, propagator_matrix=propagator_matrix, rho_old=rho_old, rho_new=rho_new, &
284                   rho_next=rho_next)
285
286      DO ispin = 1, SIZE(propagator_matrix)/2
287         re = 2*ispin - 1
288         im = 2*ispin
289
290         IF (rtp%iter == 1) THEN
291            ! For EM I have to copy rho_old onto rho_next and for ETRS,
292            ! this is the first term of the series of commutators that result in rho_next
293            CALL dbcsr_copy(rho_next(re)%matrix, rho_old(re)%matrix)
294            CALL dbcsr_copy(rho_next(im)%matrix, rho_old(im)%matrix)
295            IF (rtp_control%propagator == do_etrs) THEN
296               !since we never calculated the matrix exponential the old matrix exponential stores the unscalled propagator
297               CALL get_rtp(rtp=rtp, exp_H_old=exp_H_old, dt=dt)
298               prefac = -0.5_dp*dt
299               CALL dbcsr_scale(exp_H_old(im)%matrix, prefac)
300               IF (.NOT. rtp%do_hfx .AND. rtp_control%fixed_ions) THEN
301                  CALL bch_expansion_imaginary_propagator( &
302                     exp_H_old(im)%matrix, rho_next(re)%matrix, rho_next(im)%matrix, &
303                     rtp%filter_eps, rtp%filter_eps_small, rtp_control%eps_exp)
304               ELSE
305                  CALL dbcsr_scale(exp_H_old(re)%matrix, prefac)
306                  CALL bch_expansion_complex_propagator( &
307                     exp_H_old(re)%matrix, exp_H_old(im)%matrix, rho_next(re)%matrix, rho_next(im)%matrix, &
308                     rtp%filter_eps, rtp%filter_eps_small, rtp_control%eps_exp)
309               ENDIF
310            END IF
311         END IF
312         CALL dbcsr_copy(rho_new(re)%matrix, rho_next(re)%matrix)
313         CALL dbcsr_copy(rho_new(im)%matrix, rho_next(im)%matrix)
314         IF (.NOT. rtp%do_hfx .AND. rtp_control%fixed_ions) THEN
315            CALL bch_expansion_imaginary_propagator( &
316               propagator_matrix(im)%matrix, rho_new(re)%matrix, rho_new(im)%matrix, &
317               rtp%filter_eps, rtp%filter_eps_small, rtp_control%eps_exp)
318         ELSE
319            CALL bch_expansion_complex_propagator( &
320               propagator_matrix(re)%matrix, propagator_matrix(im)%matrix, rho_new(re)%matrix, rho_new(im)%matrix, &
321               rtp%filter_eps, rtp%filter_eps_small, rtp_control%eps_exp)
322         ENDIF
323
324      END DO
325
326      CALL timestop(handle)
327
328   END SUBROUTINE propagate_bch
329
330! **************************************************************************************************
331!> \brief decides which type of exponential has to be computed
332!> \param propagator ...
333!> \param propagator_matrix ...
334!> \param rtp_control ...
335!> \param rtp ...
336!> \author Florian Schiffmann (02.09)
337! **************************************************************************************************
338
339   SUBROUTINE compute_exponential(propagator, propagator_matrix, rtp_control, rtp)
340      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: propagator, propagator_matrix
341      TYPE(rtp_control_type), POINTER                    :: rtp_control
342      TYPE(rt_prop_type), POINTER                        :: rtp
343
344      CHARACTER(len=*), PARAMETER :: routineN = 'compute_exponential', &
345         routineP = moduleN//':'//routineN
346
347      INTEGER                                            :: i, im, ispin, re
348      TYPE(cp_fm_p_type), DIMENSION(:), POINTER          :: propagator_fm, propagator_matrix_fm
349
350      NULLIFY (propagator_fm)
351      ALLOCATE (propagator_fm(SIZE(propagator)))
352      NULLIFY (propagator_matrix_fm)
353      ALLOCATE (propagator_matrix_fm(SIZE(propagator_matrix)))
354      DO i = 1, SIZE(propagator)
355         CALL cp_fm_create(propagator_fm(i)%matrix, &
356                           matrix_struct=rtp%ao_ao_fmstruct, &
357                           name="prop_fm")
358         CALL copy_dbcsr_to_fm(propagator(i)%matrix, propagator_fm(i)%matrix)
359         CALL cp_fm_create(propagator_matrix_fm(i)%matrix, &
360                           matrix_struct=rtp%ao_ao_fmstruct, &
361                           name="prop_mat_fm")
362         CALL copy_dbcsr_to_fm(propagator_matrix(i)%matrix, propagator_matrix_fm(i)%matrix)
363      END DO
364
365      DO ispin = 1, SIZE(propagator)/2
366         re = 2*ispin - 1
367         im = 2*ispin
368
369         SELECT CASE (rtp_control%mat_exp)
370
371         CASE (do_taylor)
372            IF (rtp_control%fixed_ions .AND. .NOT. rtp%do_hfx) THEN
373               CALL taylor_only_imaginary(propagator_fm(re:im), propagator_matrix_fm(im)%matrix, &
374                                          rtp%orders(1, ispin), rtp%orders(2, ispin))
375            ELSE
376               CALL taylor_full_complex(propagator_fm(re:im), propagator_matrix_fm(re)%matrix, propagator_matrix_fm(im)%matrix, &
377                                        rtp%orders(1, ispin), rtp%orders(2, ispin))
378            END IF
379         CASE (do_pade)
380            IF (rtp_control%fixed_ions .AND. .NOT. rtp%do_hfx) THEN
381               CALL exp_pade_only_imaginary(propagator_fm(re:im), propagator_matrix_fm(im)%matrix, &
382                                            rtp%orders(1, ispin), rtp%orders(2, ispin))
383            ELSE
384               CALL exp_pade_full_complex(propagator_fm(re:im), propagator_matrix_fm(re)%matrix, propagator_matrix_fm(im)%matrix, &
385                                          rtp%orders(1, ispin), rtp%orders(2, ispin))
386            END IF
387         END SELECT
388      END DO
389
390      DO i = 1, SIZE(propagator)
391         CALL copy_fm_to_dbcsr(propagator_fm(i)%matrix, propagator(i)%matrix)
392         CALL copy_fm_to_dbcsr(propagator_matrix_fm(i)%matrix, propagator_matrix(i)%matrix)
393      END DO
394      CALL cp_fm_vect_dealloc(propagator_fm)
395      CALL cp_fm_vect_dealloc(propagator_matrix_fm)
396
397   END SUBROUTINE compute_exponential
398
399! **************************************************************************************************
400!> \brief Sparse versions of the matrix exponentials
401!> \param propagator ...
402!> \param propagator_matrix ...
403!> \param rtp_control ...
404!> \param rtp ...
405!> \author Samuel Andermatt (02.14)
406! **************************************************************************************************
407
408   SUBROUTINE compute_exponential_sparse(propagator, propagator_matrix, rtp_control, rtp)
409      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: propagator, propagator_matrix
410      TYPE(rtp_control_type), POINTER                    :: rtp_control
411      TYPE(rt_prop_type), POINTER                        :: rtp
412
413      CHARACTER(len=*), PARAMETER :: routineN = 'compute_exponential_sparse', &
414         routineP = moduleN//':'//routineN
415
416      INTEGER                                            :: handle, im, ispin, re
417
418      CALL timeset(routineN, handle)
419
420      DO ispin = 1, SIZE(propagator)/2
421         re = 2*ispin - 1
422         im = 2*ispin
423         IF (rtp_control%fixed_ions .AND. .NOT. rtp%do_hfx) THEN
424            CALL taylor_only_imaginary_dbcsr(propagator(re:im), propagator_matrix(im)%matrix, &
425                                             rtp%orders(1, ispin), rtp%orders(2, ispin), rtp%filter_eps)
426         ELSE
427            CALL taylor_full_complex_dbcsr(propagator(re:im), propagator_matrix(re)%matrix, propagator_matrix(im)%matrix, &
428                                           rtp%orders(1, ispin), rtp%orders(2, ispin), rtp%filter_eps)
429         END IF
430      END DO
431
432      CALL timestop(handle)
433
434   END SUBROUTINE compute_exponential_sparse
435
436END MODULE rt_make_propagators
437