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