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