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