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 with dbcsr matrices.
10 : !> Based on the code in matrix_exp.F from Florian Schiffmann
11 : !> \author Samuel Andermatt (02.14)
12 : ! **************************************************************************************************
13 : MODULE ls_matrix_exp
14 :
15 : USE cp_dbcsr_api, ONLY: &
16 : dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, &
17 : dbcsr_filter, dbcsr_frobenius_norm, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
18 : dbcsr_scale, dbcsr_set, dbcsr_transposed, dbcsr_type
19 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
20 : USE kinds, ONLY: dp
21 : #include "./base/base_uses.f90"
22 :
23 : IMPLICIT NONE
24 :
25 : PRIVATE
26 :
27 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ls_matrix_exp'
28 :
29 : PUBLIC :: taylor_only_imaginary_dbcsr, &
30 : taylor_full_complex_dbcsr, &
31 : cp_complex_dbcsr_gemm_3, &
32 : bch_expansion_imaginary_propagator, &
33 : bch_expansion_complex_propagator
34 :
35 : CONTAINS
36 :
37 : ! **************************************************************************************************
38 : !> \brief Convenience function. Computes the matrix multiplications needed
39 : !> for the multiplication of complex sparse matrices.
40 : !> C = beta * C + alpha * ( A ** transa ) * ( B ** transb )
41 : !> \param transa : 'N' -> normal 'T' -> transpose
42 : !> alpha,beta :: can be 0.0_dp and 1.0_dp
43 : !> \param transb ...
44 : !> \param alpha ...
45 : !> \param A_re m x k matrix ( ! for transa = 'N'), real part
46 : !> \param A_im m x k matrix ( ! for transa = 'N'), imaginary part
47 : !> \param B_re k x n matrix ( ! for transb = 'N'), real part
48 : !> \param B_im k x n matrix ( ! for transb = 'N'), imaginary part
49 : !> \param beta ...
50 : !> \param C_re m x n matrix, real part
51 : !> \param C_im m x n matrix, imaginary part
52 : !> \param filter_eps ...
53 : !> \author Samuel Andermatt
54 : !> \note
55 : !> C should have no overlap with A, B
56 : !> This subroutine uses three real matrix multiplications instead of two complex
57 : !> This reduces the amount of flops and memory bandwidth by 25%, but for memory bandwidth
58 : !> true complex algebra is still superior (one third less bandwidth needed)
59 : !> limited cases matrix multiplications
60 : ! **************************************************************************************************
61 6962 : SUBROUTINE cp_complex_dbcsr_gemm_3(transa, transb, alpha, A_re, A_im, &
62 : B_re, B_im, beta, C_re, C_im, filter_eps)
63 : CHARACTER(LEN=1), INTENT(IN) :: transa, transb
64 : REAL(KIND=dp), INTENT(IN) :: alpha
65 : TYPE(dbcsr_type), INTENT(IN) :: A_re, A_im, B_re, B_im
66 : REAL(KIND=dp), INTENT(IN) :: beta
67 : TYPE(dbcsr_type), INTENT(INOUT) :: C_re, C_im
68 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: filter_eps
69 :
70 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_complex_dbcsr_gemm_3'
71 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
72 :
73 : CHARACTER(LEN=1) :: transa2, transb2
74 : INTEGER :: handle
75 : REAL(KIND=dp) :: alpha2, alpha3, alpha4
76 : TYPE(dbcsr_type), POINTER :: a_plus_b, ac, bd, c_plus_d
77 :
78 6962 : CALL timeset(routineN, handle)
79 : !A complex matrix matrix multiplication can be done with only three multiplications
80 : !(a+ib)*(c+id)=ac-bd+i((a+b)*(c+d) - ac - bd)
81 : !A_re=a, A_im=b, B_re=c, B_im=d
82 :
83 6962 : alpha2 = -alpha
84 6962 : alpha3 = alpha
85 6962 : alpha4 = alpha
86 :
87 6962 : IF (transa == "C") THEN
88 0 : alpha2 = -alpha2
89 0 : alpha3 = -alpha3
90 0 : transa2 = "T"
91 : ELSE
92 6962 : transa2 = transa
93 : END IF
94 6962 : IF (transb == "C") THEN
95 1072 : alpha2 = -alpha2
96 1072 : alpha4 = -alpha4
97 1072 : transb2 = "T"
98 : ELSE
99 5890 : transb2 = transb
100 : END IF
101 :
102 : !create the work matrices
103 : NULLIFY (ac)
104 6962 : ALLOCATE (ac)
105 6962 : CALL dbcsr_create(ac, template=A_re, matrix_type="N")
106 : NULLIFY (bd)
107 6962 : ALLOCATE (bd)
108 6962 : CALL dbcsr_create(bd, template=A_re, matrix_type="N")
109 : NULLIFY (a_plus_b)
110 6962 : ALLOCATE (a_plus_b)
111 6962 : CALL dbcsr_create(a_plus_b, template=A_re, matrix_type="N")
112 : NULLIFY (c_plus_d)
113 6962 : ALLOCATE (c_plus_d)
114 6962 : CALL dbcsr_create(c_plus_d, template=A_re, matrix_type="N")
115 :
116 : !Do the neccesarry multiplications
117 6962 : CALL dbcsr_multiply(transa2, transb2, alpha, A_re, B_re, zero, ac, filter_eps=filter_eps)
118 6962 : CALL dbcsr_multiply(transa2, transb2, alpha2, A_im, B_im, zero, bd, filter_eps=filter_eps)
119 :
120 6962 : CALL dbcsr_add(a_plus_b, A_re, zero, alpha)
121 6962 : CALL dbcsr_add(a_plus_b, A_im, one, alpha3)
122 6962 : CALL dbcsr_add(c_plus_d, B_re, zero, alpha)
123 6962 : CALL dbcsr_add(c_plus_d, B_im, one, alpha4)
124 :
125 : !this can already be written into C_im
126 : !now both matrixes have been scaled which means we currently multiplied by alpha squared
127 6962 : CALL dbcsr_multiply(transa2, transb2, one/alpha, a_plus_b, c_plus_d, beta, C_im, filter_eps=filter_eps)
128 :
129 : !now add up all the terms into the result
130 6962 : CALL dbcsr_add(C_re, ac, beta, one)
131 : !the minus sign was already taken care of at the definition of alpha2
132 6962 : CALL dbcsr_add(C_re, bd, one, one)
133 6962 : CALL dbcsr_filter(C_re, filter_eps)
134 :
135 6962 : CALL dbcsr_add(C_im, ac, one, -one)
136 : !the minus sign was already taken care of at the definition of alpha2
137 6962 : CALL dbcsr_add(C_im, bd, one, one)
138 6962 : CALL dbcsr_filter(C_im, filter_eps)
139 :
140 : !Deallocate the work matrices
141 6962 : CALL dbcsr_deallocate_matrix(ac)
142 6962 : CALL dbcsr_deallocate_matrix(bd)
143 6962 : CALL dbcsr_deallocate_matrix(a_plus_b)
144 6962 : CALL dbcsr_deallocate_matrix(c_plus_d)
145 :
146 6962 : CALL timestop(handle)
147 :
148 6962 : END SUBROUTINE cp_complex_dbcsr_gemm_3
149 :
150 : ! **************************************************************************************************
151 : !> \brief specialized subroutine for purely imaginary matrix exponentials
152 : !> \param exp_H ...
153 : !> \param im_matrix ...
154 : !> \param nsquare ...
155 : !> \param ntaylor ...
156 : !> \param filter_eps ...
157 : !> \author Samuel Andermatt (02.2014)
158 : ! **************************************************************************************************
159 708 : SUBROUTINE taylor_only_imaginary_dbcsr(exp_H, im_matrix, nsquare, ntaylor, filter_eps)
160 :
161 : TYPE(dbcsr_p_type), DIMENSION(2) :: exp_H
162 : TYPE(dbcsr_type), POINTER :: im_matrix
163 : INTEGER, INTENT(in) :: nsquare, ntaylor
164 : REAL(KIND=dp), INTENT(in) :: filter_eps
165 :
166 : CHARACTER(len=*), PARAMETER :: routineN = 'taylor_only_imaginary_dbcsr'
167 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
168 :
169 : INTEGER :: handle, i, nloop
170 : REAL(KIND=dp) :: square_fac, Tfac, tmp
171 : TYPE(dbcsr_type), POINTER :: T1, T2, Tres_im, Tres_re
172 :
173 708 : CALL timeset(routineN, handle)
174 :
175 : !The divider that comes from the scaling and squaring procedure
176 708 : square_fac = 1.0_dp/(2.0_dp**REAL(nsquare, dp))
177 :
178 : !Allocate work matrices
179 : NULLIFY (T1)
180 708 : ALLOCATE (T1)
181 708 : CALL dbcsr_create(T1, template=im_matrix, matrix_type="N")
182 : NULLIFY (T2)
183 708 : ALLOCATE (T2)
184 708 : CALL dbcsr_create(T2, template=im_matrix, matrix_type="N")
185 : NULLIFY (Tres_re)
186 708 : ALLOCATE (Tres_re)
187 708 : CALL dbcsr_create(Tres_re, template=im_matrix, matrix_type="N")
188 : NULLIFY (Tres_im)
189 708 : ALLOCATE (Tres_im)
190 708 : CALL dbcsr_create(Tres_im, template=im_matrix, matrix_type="N")
191 :
192 : !Create the unit matrices
193 708 : CALL dbcsr_set(T1, zero)
194 708 : CALL dbcsr_add_on_diag(T1, one)
195 708 : CALL dbcsr_set(Tres_re, zero)
196 708 : CALL dbcsr_add_on_diag(Tres_re, one)
197 708 : CALL dbcsr_set(Tres_im, zero)
198 :
199 708 : nloop = CEILING(REAL(ntaylor, dp)/2.0_dp)
200 : !the inverse of the prefactor in the taylor series
201 708 : tmp = 1.0_dp
202 3516 : DO i = 1, nloop
203 2808 : CALL dbcsr_scale(T1, 1.0_dp/(REAL(i, dp)*2.0_dp - 1.0_dp))
204 2808 : CALL dbcsr_filter(T1, filter_eps)
205 : CALL dbcsr_multiply("N", "N", square_fac, im_matrix, T1, zero, &
206 2808 : T2, filter_eps=filter_eps)
207 2808 : Tfac = one
208 2808 : IF (MOD(i, 2) == 0) Tfac = -Tfac
209 2808 : CALL dbcsr_add(Tres_im, T2, one, Tfac)
210 2808 : CALL dbcsr_scale(T2, 1.0_dp/(REAL(i, dp)*2.0_dp))
211 2808 : CALL dbcsr_filter(T2, filter_eps)
212 : CALL dbcsr_multiply("N", "N", square_fac, im_matrix, T2, zero, &
213 2808 : T1, filter_eps=filter_eps)
214 2808 : Tfac = one
215 2808 : IF (MOD(i, 2) == 1) Tfac = -Tfac
216 3516 : CALL dbcsr_add(Tres_re, T1, one, Tfac)
217 : END DO
218 :
219 : !Square the matrices, due to the scaling and squaring procedure
220 708 : IF (nsquare .GT. 0) THEN
221 0 : DO i = 1, nsquare
222 : CALL cp_complex_dbcsr_gemm_3("N", "N", one, Tres_re, Tres_im, &
223 : Tres_re, Tres_im, zero, exp_H(1)%matrix, exp_H(2)%matrix, &
224 0 : filter_eps=filter_eps)
225 0 : CALL dbcsr_copy(Tres_re, exp_H(1)%matrix)
226 0 : CALL dbcsr_copy(Tres_im, exp_H(2)%matrix)
227 : END DO
228 : ELSE
229 708 : CALL dbcsr_copy(exp_H(1)%matrix, Tres_re)
230 708 : CALL dbcsr_copy(exp_H(2)%matrix, Tres_im)
231 : END IF
232 708 : CALL dbcsr_deallocate_matrix(T1)
233 708 : CALL dbcsr_deallocate_matrix(T2)
234 708 : CALL dbcsr_deallocate_matrix(Tres_re)
235 708 : CALL dbcsr_deallocate_matrix(Tres_im)
236 :
237 708 : CALL timestop(handle)
238 :
239 708 : END SUBROUTINE taylor_only_imaginary_dbcsr
240 :
241 : ! **************************************************************************************************
242 : !> \brief subroutine for general complex matrix exponentials
243 : !> on input a separate dbcsr_type for real and complex part
244 : !> on output a size 2 dbcsr_p_type, first element is the real part of
245 : !> the exponential second the imaginary
246 : !> \param exp_H ...
247 : !> \param re_part ...
248 : !> \param im_part ...
249 : !> \param nsquare ...
250 : !> \param ntaylor ...
251 : !> \param filter_eps ...
252 : !> \author Samuel Andermatt (02.2014)
253 : ! **************************************************************************************************
254 264 : SUBROUTINE taylor_full_complex_dbcsr(exp_H, re_part, im_part, nsquare, ntaylor, filter_eps)
255 : TYPE(dbcsr_p_type), DIMENSION(2) :: exp_H
256 : TYPE(dbcsr_type), POINTER :: re_part, im_part
257 : INTEGER, INTENT(in) :: nsquare, ntaylor
258 : REAL(KIND=dp), INTENT(in) :: filter_eps
259 :
260 : CHARACTER(len=*), PARAMETER :: routineN = 'taylor_full_complex_dbcsr'
261 :
262 : INTEGER :: handle, i
263 : REAL(KIND=dp) :: square_fac
264 : TYPE(dbcsr_type) :: T1_im, T1_re, T2_im, T2_re
265 : TYPE(dbcsr_type), POINTER :: Tres_im, Tres_re
266 :
267 264 : CALL timeset(routineN, handle)
268 :
269 : ! convenient aliases for result matrices
270 264 : Tres_re => exp_H(1)%matrix
271 264 : Tres_im => exp_H(2)%matrix
272 :
273 : ! The divider that comes from the scaling and squaring procedure
274 264 : square_fac = 1.0_dp/REAL(2**nsquare, dp)
275 :
276 : ! Allocate work matrices
277 264 : CALL dbcsr_create(T1_re, template=re_part, matrix_type="N")
278 264 : CALL dbcsr_create(T1_im, template=re_part, matrix_type="N")
279 264 : CALL dbcsr_create(T2_re, template=re_part, matrix_type="N")
280 264 : CALL dbcsr_create(T2_im, template=re_part, matrix_type="N")
281 :
282 : ! T1 = identity
283 264 : CALL dbcsr_set(T1_re, 0.0_dp)
284 264 : CALL dbcsr_set(T1_im, 0.0_dp)
285 264 : CALL dbcsr_add_on_diag(T1_re, 1.0_dp)
286 :
287 : ! Tres = identity
288 264 : CALL dbcsr_set(Tres_re, 0.0_dp)
289 264 : CALL dbcsr_set(Tres_im, 0.0_dp)
290 264 : CALL dbcsr_add_on_diag(Tres_re, 1.0_dp)
291 :
292 2522 : DO i = 1, ntaylor
293 : ! T1 = T1 / i
294 2258 : CALL dbcsr_scale(T1_re, 1.0_dp/REAL(i, dp))
295 2258 : CALL dbcsr_scale(T1_im, 1.0_dp/REAL(i, dp))
296 2258 : CALL dbcsr_filter(T1_re, filter_eps)
297 2258 : CALL dbcsr_filter(T1_im, filter_eps)
298 :
299 : ! T2 = square_fac * T1 * input
300 : CALL cp_complex_dbcsr_gemm_3("N", "N", alpha=square_fac, beta=0.0_dp, &
301 : A_re=T1_re, A_im=T1_im, &
302 : B_re=re_part, B_im=im_part, &
303 2258 : C_re=T2_re, C_im=T2_im, filter_eps=filter_eps)
304 :
305 : ! Tres = Tres + T2
306 2258 : CALL dbcsr_add(Tres_re, T2_re, 1.0_dp, 1.0_dp)
307 2258 : CALL dbcsr_add(Tres_im, T2_im, 1.0_dp, 1.0_dp)
308 :
309 : ! T1 = T2
310 2258 : CALL dbcsr_copy(T1_re, T2_re)
311 2522 : CALL dbcsr_copy(T1_im, T2_im)
312 : END DO
313 :
314 264 : IF (nsquare .GT. 0) THEN
315 970 : DO i = 1, nsquare
316 : ! T2 = Tres * Tres
317 : CALL cp_complex_dbcsr_gemm_3("N", "N", alpha=1.0_dp, beta=0.0_dp, &
318 : A_re=Tres_re, A_im=Tres_im, &
319 : B_re=Tres_re, B_im=Tres_im, &
320 758 : C_re=T2_re, C_im=T2_im, filter_eps=filter_eps)
321 :
322 : ! Tres = T2
323 758 : CALL dbcsr_copy(Tres_re, T2_re)
324 970 : CALL dbcsr_copy(Tres_im, T2_im)
325 : END DO
326 : END IF
327 :
328 264 : CALL dbcsr_release(T1_re)
329 264 : CALL dbcsr_release(T1_im)
330 264 : CALL dbcsr_release(T2_re)
331 264 : CALL dbcsr_release(T2_im)
332 :
333 264 : CALL timestop(handle)
334 :
335 264 : END SUBROUTINE taylor_full_complex_dbcsr
336 :
337 : ! **************************************************************************************************
338 : !> \brief The Baker-Campbell-Hausdorff expansion for a purely imaginary exponent (e.g. rtp)
339 : !> Works for a non unitary propagator, because the density matrix is hermitian
340 : !> \param propagator The exponent of the matrix exponential
341 : !> \param density_re Real part of the density matrix
342 : !> \param density_im Imaginary part of the density matrix
343 : !> \param filter_eps The filtering threshold for all matrices
344 : !> \param filter_eps_small ...
345 : !> \param eps_exp The accuracy of the exponential
346 : !> \author Samuel Andermatt (02.2014)
347 : ! **************************************************************************************************
348 112 : SUBROUTINE bch_expansion_imaginary_propagator(propagator, density_re, density_im, filter_eps, filter_eps_small, eps_exp)
349 : TYPE(dbcsr_type), POINTER :: propagator, density_re, density_im
350 : REAL(KIND=dp), INTENT(in) :: filter_eps, filter_eps_small, eps_exp
351 :
352 : CHARACTER(len=*), PARAMETER :: routineN = 'bch_expansion_imaginary_propagator'
353 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
354 :
355 : INTEGER :: handle, i, unit_nr
356 : LOGICAL :: convergence
357 : REAL(KIND=dp) :: alpha, max_alpha, prefac
358 : TYPE(dbcsr_type), POINTER :: comm, comm2, tmp, tmp2
359 :
360 112 : CALL timeset(routineN, handle)
361 :
362 112 : unit_nr = cp_logger_get_default_io_unit()
363 :
364 : NULLIFY (tmp)
365 112 : ALLOCATE (tmp)
366 112 : CALL dbcsr_create(tmp, template=propagator)
367 : NULLIFY (tmp2)
368 112 : ALLOCATE (tmp2)
369 112 : CALL dbcsr_create(tmp2, template=propagator)
370 : NULLIFY (comm)
371 112 : ALLOCATE (comm)
372 112 : CALL dbcsr_create(comm, template=propagator)
373 : NULLIFY (comm2)
374 112 : ALLOCATE (comm2)
375 112 : CALL dbcsr_create(comm2, template=propagator)
376 :
377 112 : CALL dbcsr_copy(tmp, density_re)
378 112 : CALL dbcsr_copy(tmp2, density_im)
379 :
380 112 : convergence = .FALSE.
381 494 : DO i = 1, 20
382 494 : prefac = one/i
383 : CALL dbcsr_multiply("N", "N", -prefac, propagator, tmp2, zero, comm, &
384 494 : filter_eps=filter_eps_small)
385 : CALL dbcsr_multiply("N", "N", prefac, propagator, tmp, zero, comm2, &
386 494 : filter_eps=filter_eps_small)
387 494 : CALL dbcsr_transposed(tmp, comm)
388 494 : CALL dbcsr_transposed(tmp2, comm2)
389 494 : CALL dbcsr_add(comm, tmp, one, one)
390 494 : CALL dbcsr_add(comm2, tmp2, one, -one)
391 494 : CALL dbcsr_add(density_re, comm, one, one)
392 494 : CALL dbcsr_add(density_im, comm2, one, one)
393 494 : CALL dbcsr_copy(tmp, comm)
394 494 : CALL dbcsr_copy(tmp2, comm2)
395 : !check for convergence
396 494 : max_alpha = zero
397 494 : alpha = dbcsr_frobenius_norm(comm)
398 494 : IF (alpha > max_alpha) max_alpha = alpha
399 494 : alpha = dbcsr_frobenius_norm(comm2)
400 494 : IF (alpha > max_alpha) max_alpha = alpha
401 494 : IF (max_alpha < eps_exp) convergence = .TRUE.
402 382 : IF (convergence) THEN
403 112 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="((T3,A,I2,A))") &
404 56 : "BCH converged after ", i, " steps"
405 : EXIT
406 : END IF
407 : END DO
408 :
409 112 : CALL dbcsr_filter(density_re, filter_eps)
410 112 : CALL dbcsr_filter(density_im, filter_eps)
411 :
412 112 : CPWARN_IF(.NOT. convergence, "BCH method did not converge")
413 :
414 112 : CALL dbcsr_deallocate_matrix(tmp)
415 112 : CALL dbcsr_deallocate_matrix(tmp2)
416 112 : CALL dbcsr_deallocate_matrix(comm)
417 112 : CALL dbcsr_deallocate_matrix(comm2)
418 :
419 112 : CALL timestop(handle)
420 :
421 112 : END SUBROUTINE bch_expansion_imaginary_propagator
422 :
423 : ! **************************************************************************************************
424 : !> \brief The Baker-Campbell-Hausdorff expansion for a complex exponent (e.g. rtp)
425 : !> Works for a non unitary propagator, because the density matrix is hermitian
426 : !> \param propagator_re Real part of the exponent
427 : !> \param propagator_im Imaginary part of the exponent
428 : !> \param density_re Real part of the density matrix
429 : !> \param density_im Imaginary part of the density matrix
430 : !> \param filter_eps The filtering threshold for all matrices
431 : !> \param filter_eps_small ...
432 : !> \param eps_exp The accuracy of the exponential
433 : !> \author Samuel Andermatt (02.2014)
434 : ! **************************************************************************************************
435 130 : SUBROUTINE bch_expansion_complex_propagator(propagator_re, propagator_im, density_re, density_im, &
436 : filter_eps, filter_eps_small, eps_exp)
437 : TYPE(dbcsr_type), POINTER :: propagator_re, propagator_im, &
438 : density_re, density_im
439 : REAL(KIND=dp), INTENT(in) :: filter_eps, filter_eps_small, eps_exp
440 :
441 : CHARACTER(len=*), PARAMETER :: routineN = 'bch_expansion_complex_propagator'
442 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
443 :
444 : INTEGER :: handle, i, unit_nr
445 : LOGICAL :: convergence
446 : REAL(KIND=dp) :: alpha, max_alpha, prefac
447 : TYPE(dbcsr_type), POINTER :: comm, comm2, tmp, tmp2
448 :
449 130 : CALL timeset(routineN, handle)
450 :
451 130 : unit_nr = cp_logger_get_default_io_unit()
452 :
453 : NULLIFY (tmp)
454 130 : ALLOCATE (tmp)
455 130 : CALL dbcsr_create(tmp, template=propagator_re)
456 : NULLIFY (tmp2)
457 130 : ALLOCATE (tmp2)
458 130 : CALL dbcsr_create(tmp2, template=propagator_re)
459 : NULLIFY (comm)
460 130 : ALLOCATE (comm)
461 130 : CALL dbcsr_create(comm, template=propagator_re)
462 : NULLIFY (comm2)
463 130 : ALLOCATE (comm2)
464 130 : CALL dbcsr_create(comm2, template=propagator_re)
465 :
466 130 : CALL dbcsr_copy(tmp, density_re)
467 130 : CALL dbcsr_copy(tmp2, density_im)
468 :
469 130 : convergence = .FALSE.
470 :
471 1254 : DO i = 1, 20
472 1254 : prefac = one/i
473 : CALL cp_complex_dbcsr_gemm_3("N", "N", prefac, propagator_re, propagator_im, &
474 1254 : tmp, tmp2, zero, comm, comm2, filter_eps=filter_eps_small)
475 1254 : CALL dbcsr_transposed(tmp, comm)
476 1254 : CALL dbcsr_transposed(tmp2, comm2)
477 1254 : CALL dbcsr_add(comm, tmp, one, one)
478 1254 : CALL dbcsr_add(comm2, tmp2, one, -one)
479 1254 : CALL dbcsr_add(density_re, comm, one, one)
480 1254 : CALL dbcsr_add(density_im, comm2, one, one)
481 1254 : CALL dbcsr_copy(tmp, comm)
482 1254 : CALL dbcsr_copy(tmp2, comm2)
483 : !check for convergence
484 1254 : max_alpha = zero
485 1254 : alpha = dbcsr_frobenius_norm(comm)
486 1254 : IF (alpha > max_alpha) max_alpha = alpha
487 1254 : alpha = dbcsr_frobenius_norm(comm2)
488 1254 : IF (alpha > max_alpha) max_alpha = alpha
489 1254 : IF (max_alpha < eps_exp) convergence = .TRUE.
490 1124 : IF (convergence) THEN
491 130 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="((T3,A,I2,A))") &
492 65 : "BCH converged after ", i, " steps"
493 : EXIT
494 : END IF
495 : END DO
496 :
497 130 : CALL dbcsr_filter(density_re, filter_eps)
498 130 : CALL dbcsr_filter(density_im, filter_eps)
499 :
500 130 : CPWARN_IF(.NOT. convergence, "BCH method did not converge")
501 :
502 130 : CALL dbcsr_deallocate_matrix(tmp)
503 130 : CALL dbcsr_deallocate_matrix(tmp2)
504 130 : CALL dbcsr_deallocate_matrix(comm)
505 130 : CALL dbcsr_deallocate_matrix(comm2)
506 :
507 130 : CALL timestop(handle)
508 :
509 130 : END SUBROUTINE bch_expansion_complex_propagator
510 :
511 : END MODULE ls_matrix_exp
|
