Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 :
14 : MODULE ls_matrix_exp
15 :
16 : USE cp_dbcsr_api, ONLY: &
17 : dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, &
18 : dbcsr_filter, dbcsr_frobenius_norm, dbcsr_multiply, dbcsr_p_type, dbcsr_scale, dbcsr_set, &
19 : dbcsr_transposed, dbcsr_type, dbcsr_type_complex_8
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 :
63 3946 : 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 3946 : 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 3946 : alpha2 = -alpha
86 3946 : alpha3 = alpha
87 3946 : alpha4 = alpha
88 :
89 3946 : IF (transa == "C") THEN
90 0 : alpha2 = -alpha2
91 0 : alpha3 = -alpha3
92 0 : transa2 = "T"
93 : ELSE
94 3946 : transa2 = transa
95 : END IF
96 3946 : IF (transb == "C") THEN
97 1072 : alpha2 = -alpha2
98 1072 : alpha4 = -alpha4
99 1072 : transb2 = "T"
100 : ELSE
101 2874 : transb2 = transb
102 : END IF
103 :
104 : !create the work matrices
105 : NULLIFY (ac)
106 3946 : ALLOCATE (ac)
107 3946 : CALL dbcsr_create(ac, template=A_re, matrix_type="N")
108 : NULLIFY (bd)
109 3946 : ALLOCATE (bd)
110 3946 : CALL dbcsr_create(bd, template=A_re, matrix_type="N")
111 : NULLIFY (a_plus_b)
112 3946 : ALLOCATE (a_plus_b)
113 3946 : CALL dbcsr_create(a_plus_b, template=A_re, matrix_type="N")
114 : NULLIFY (c_plus_d)
115 3946 : ALLOCATE (c_plus_d)
116 3946 : CALL dbcsr_create(c_plus_d, template=A_re, matrix_type="N")
117 :
118 : !Do the neccesarry multiplications
119 3946 : CALL dbcsr_multiply(transa2, transb2, alpha, A_re, B_re, zero, ac, filter_eps=filter_eps)
120 3946 : CALL dbcsr_multiply(transa2, transb2, alpha2, A_im, B_im, zero, bd, filter_eps=filter_eps)
121 :
122 3946 : CALL dbcsr_add(a_plus_b, A_re, zero, alpha)
123 3946 : CALL dbcsr_add(a_plus_b, A_im, one, alpha3)
124 3946 : CALL dbcsr_add(c_plus_d, B_re, zero, alpha)
125 3946 : 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 3946 : 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 3946 : CALL dbcsr_add(C_re, ac, beta, one)
133 : !the minus sign was already taken care of at the definition of alpha2
134 3946 : CALL dbcsr_add(C_re, bd, one, one)
135 3946 : CALL dbcsr_filter(C_re, filter_eps)
136 :
137 3946 : CALL dbcsr_add(C_im, ac, one, -one)
138 : !the minus sign was already taken care of at the definition of alpha2
139 3946 : CALL dbcsr_add(C_im, bd, one, one)
140 3946 : CALL dbcsr_filter(C_im, filter_eps)
141 :
142 : !Deallocate the work matrices
143 3946 : CALL dbcsr_deallocate_matrix(ac)
144 3946 : CALL dbcsr_deallocate_matrix(bd)
145 3946 : CALL dbcsr_deallocate_matrix(a_plus_b)
146 3946 : CALL dbcsr_deallocate_matrix(c_plus_d)
147 :
148 3946 : CALL timestop(handle)
149 :
150 3946 : END SUBROUTINE
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 :
162 708 : SUBROUTINE taylor_only_imaginary_dbcsr(exp_H, im_matrix, nsquare, ntaylor, filter_eps)
163 :
164 : TYPE(dbcsr_p_type), DIMENSION(2) :: exp_H
165 : TYPE(dbcsr_type), POINTER :: im_matrix
166 : INTEGER, INTENT(in) :: nsquare, ntaylor
167 : REAL(KIND=dp), INTENT(in) :: filter_eps
168 :
169 : CHARACTER(len=*), PARAMETER :: routineN = 'taylor_only_imaginary_dbcsr'
170 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
171 :
172 : INTEGER :: handle, i, nloop
173 : REAL(KIND=dp) :: square_fac, Tfac, tmp
174 : TYPE(dbcsr_type), POINTER :: T1, T2, Tres_im, Tres_re
175 :
176 708 : CALL timeset(routineN, handle)
177 :
178 : !The divider that comes from the scaling and squaring procedure
179 708 : square_fac = 1.0_dp/(2.0_dp**REAL(nsquare, dp))
180 :
181 : !Allocate work matrices
182 : NULLIFY (T1)
183 708 : ALLOCATE (T1)
184 708 : CALL dbcsr_create(T1, template=im_matrix, matrix_type="N")
185 : NULLIFY (T2)
186 708 : ALLOCATE (T2)
187 708 : CALL dbcsr_create(T2, template=im_matrix, matrix_type="N")
188 : NULLIFY (Tres_re)
189 708 : ALLOCATE (Tres_re)
190 708 : CALL dbcsr_create(Tres_re, template=im_matrix, matrix_type="N")
191 : NULLIFY (Tres_im)
192 708 : ALLOCATE (Tres_im)
193 708 : CALL dbcsr_create(Tres_im, template=im_matrix, matrix_type="N")
194 :
195 : !Create the unit matrices
196 708 : CALL dbcsr_set(T1, zero)
197 708 : CALL dbcsr_add_on_diag(T1, one)
198 708 : CALL dbcsr_set(Tres_re, zero)
199 708 : CALL dbcsr_add_on_diag(Tres_re, one)
200 708 : CALL dbcsr_set(Tres_im, zero)
201 :
202 708 : nloop = CEILING(REAL(ntaylor, dp)/2.0_dp)
203 : !the inverse of the prefactor in the taylor series
204 708 : tmp = 1.0_dp
205 3516 : DO i = 1, nloop
206 2808 : CALL dbcsr_scale(T1, 1.0_dp/(REAL(i, dp)*2.0_dp - 1.0_dp))
207 2808 : CALL dbcsr_filter(T1, filter_eps)
208 : CALL dbcsr_multiply("N", "N", square_fac, im_matrix, T1, zero, &
209 2808 : T2, filter_eps=filter_eps)
210 2808 : Tfac = one
211 2808 : IF (MOD(i, 2) == 0) Tfac = -Tfac
212 2808 : CALL dbcsr_add(Tres_im, T2, one, Tfac)
213 2808 : CALL dbcsr_scale(T2, 1.0_dp/(REAL(i, dp)*2.0_dp))
214 2808 : CALL dbcsr_filter(T2, filter_eps)
215 : CALL dbcsr_multiply("N", "N", square_fac, im_matrix, T2, zero, &
216 2808 : T1, filter_eps=filter_eps)
217 2808 : Tfac = one
218 2808 : IF (MOD(i, 2) == 1) Tfac = -Tfac
219 3516 : CALL dbcsr_add(Tres_re, T1, one, Tfac)
220 : END DO
221 :
222 : !Square the matrices, due to the scaling and squaring procedure
223 708 : IF (nsquare .GT. 0) THEN
224 0 : DO i = 1, nsquare
225 : CALL cp_complex_dbcsr_gemm_3("N", "N", one, Tres_re, Tres_im, &
226 : Tres_re, Tres_im, zero, exp_H(1)%matrix, exp_H(2)%matrix, &
227 0 : filter_eps=filter_eps)
228 0 : CALL dbcsr_copy(Tres_re, exp_H(1)%matrix)
229 0 : CALL dbcsr_copy(Tres_im, exp_H(2)%matrix)
230 : END DO
231 : ELSE
232 708 : CALL dbcsr_copy(exp_H(1)%matrix, Tres_re)
233 708 : CALL dbcsr_copy(exp_H(2)%matrix, Tres_im)
234 : END IF
235 708 : CALL dbcsr_deallocate_matrix(T1)
236 708 : CALL dbcsr_deallocate_matrix(T2)
237 708 : CALL dbcsr_deallocate_matrix(Tres_re)
238 708 : CALL dbcsr_deallocate_matrix(Tres_im)
239 :
240 708 : CALL timestop(handle)
241 :
242 708 : END SUBROUTINE taylor_only_imaginary_dbcsr
243 :
244 : ! **************************************************************************************************
245 : !> \brief subroutine for general complex matrix exponentials
246 : !> on input a separate dbcsr_type for real and complex part
247 : !> on output a size 2 dbcsr_p_type, first element is the real part of
248 : !> the exponential second the imaginary
249 : !> \param exp_H ...
250 : !> \param re_part ...
251 : !> \param im_part ...
252 : !> \param nsquare ...
253 : !> \param ntaylor ...
254 : !> \param filter_eps ...
255 : !> \author Samuel Andermatt (02.2014)
256 : ! **************************************************************************************************
257 264 : SUBROUTINE taylor_full_complex_dbcsr(exp_H, re_part, im_part, nsquare, ntaylor, filter_eps)
258 : TYPE(dbcsr_p_type), DIMENSION(2) :: exp_H
259 : TYPE(dbcsr_type), POINTER :: re_part, im_part
260 : INTEGER, INTENT(in) :: nsquare, ntaylor
261 : REAL(KIND=dp), INTENT(in) :: filter_eps
262 :
263 : CHARACTER(len=*), PARAMETER :: routineN = 'taylor_full_complex_dbcsr'
264 : COMPLEX(KIND=dp), PARAMETER :: one = (1.0_dp, 0.0_dp), &
265 : zero = (0.0_dp, 0.0_dp)
266 :
267 : COMPLEX(KIND=dp) :: square_fac
268 : INTEGER :: handle, i
269 : TYPE(dbcsr_type), POINTER :: T1, T2, T3, Tres
270 :
271 264 : CALL timeset(routineN, handle)
272 :
273 : !The divider that comes from the scaling and squaring procedure
274 264 : square_fac = CMPLX(1.0_dp/(2.0_dp**REAL(nsquare, dp)), 0.0_dp, KIND=dp)
275 :
276 : !Allocate work matrices
277 : NULLIFY (T1)
278 264 : ALLOCATE (T1)
279 : CALL dbcsr_create(T1, template=re_part, matrix_type="N", &
280 264 : data_type=dbcsr_type_complex_8)
281 : NULLIFY (T2)
282 264 : ALLOCATE (T2)
283 : CALL dbcsr_create(T2, template=re_part, matrix_type="N", &
284 264 : data_type=dbcsr_type_complex_8)
285 : NULLIFY (T3)
286 264 : ALLOCATE (T3)
287 : CALL dbcsr_create(T3, template=re_part, matrix_type="N", &
288 264 : data_type=dbcsr_type_complex_8)
289 : NULLIFY (Tres)
290 264 : ALLOCATE (Tres)
291 : CALL dbcsr_create(Tres, template=re_part, matrix_type="N", &
292 264 : data_type=dbcsr_type_complex_8)
293 :
294 : !Fuse the input matrices to a single complex matrix
295 264 : CALL dbcsr_copy(T3, re_part)
296 264 : CALL dbcsr_copy(Tres, im_part) !will later on be set back to zero
297 264 : CALL dbcsr_scale(Tres, CMPLX(0.0_dp, 1.0_dp, KIND=dp))
298 264 : CALL dbcsr_add(T3, Tres, one, one)
299 :
300 : !Create the unit matrices
301 264 : CALL dbcsr_set(T1, zero)
302 264 : CALL dbcsr_add_on_diag(T1, one)
303 264 : CALL dbcsr_set(Tres, zero)
304 264 : CALL dbcsr_add_on_diag(Tres, one)
305 :
306 2522 : DO i = 1, ntaylor
307 2258 : CALL dbcsr_scale(T1, one/CMPLX(i*1.0_dp, 0.0_dp, KIND=dp))
308 2258 : CALL dbcsr_filter(T1, filter_eps)
309 : CALL dbcsr_multiply("N", "N", square_fac, T1, T3, &
310 2258 : zero, T2, filter_eps=filter_eps)
311 2258 : CALL dbcsr_add(Tres, T2, one, one)
312 2522 : CALL dbcsr_copy(T1, T2)
313 : END DO
314 :
315 264 : IF (nsquare .GT. 0) THEN
316 970 : DO i = 1, nsquare
317 : CALL dbcsr_multiply("N", "N", one, Tres, Tres, zero, &
318 758 : T2, filter_eps=filter_eps)
319 970 : CALL dbcsr_copy(Tres, T2)
320 : END DO
321 : END IF
322 :
323 264 : CALL dbcsr_copy(exp_H(1)%matrix, Tres, keep_imaginary=.FALSE.)
324 264 : CALL dbcsr_scale(Tres, CMPLX(0.0_dp, -1.0_dp, KIND=dp))
325 264 : CALL dbcsr_copy(exp_H(2)%matrix, Tres, keep_imaginary=.FALSE.)
326 :
327 264 : CALL dbcsr_deallocate_matrix(T1)
328 264 : CALL dbcsr_deallocate_matrix(T2)
329 264 : CALL dbcsr_deallocate_matrix(T3)
330 264 : CALL dbcsr_deallocate_matrix(Tres)
331 :
332 264 : CALL timestop(handle)
333 :
334 264 : END SUBROUTINE taylor_full_complex_dbcsr
335 :
336 : ! **************************************************************************************************
337 : !> \brief The Baker-Campbell-Hausdorff expansion for a purely imaginary exponent (e.g. rtp)
338 : !> Works for a non unitary propagator, because the density matrix is hermitian
339 : !> \param propagator The exponent of the matrix exponential
340 : !> \param density_re Real part of the density matrix
341 : !> \param density_im Imaginary part of the density matrix
342 : !> \param filter_eps The filtering threshold for all matrices
343 : !> \param filter_eps_small ...
344 : !> \param eps_exp The accuracy of the exponential
345 : !> \author Samuel Andermatt (02.2014)
346 : ! **************************************************************************************************
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
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 :
436 130 : SUBROUTINE bch_expansion_complex_propagator(propagator_re, propagator_im, density_re, density_im, filter_eps, &
437 : 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
511 :
512 : END MODULE ls_matrix_exp
|
