(git:936074a)
Loading...
Searching...
No Matches
rt_propagator_init.F
Go to the documentation of this file.
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 that prepare rtp and EMD
10!> \author Florian Schiffmann (02.09)
11! **************************************************************************************************
12MODULE rt_propagator_init
13
14 USE arnoldi_api, ONLY: arnoldi_extremal
15 USE cp_control_types, ONLY: dft_control_type,&
16 rtp_control_type
17 USE cp_dbcsr_api, ONLY: &
18 dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_filter, dbcsr_multiply, &
19 dbcsr_p_type, dbcsr_scale, dbcsr_set, dbcsr_type
20 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
21 copy_fm_to_dbcsr,&
22 cp_dbcsr_plus_fm_fm_t
23 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
24 cp_fm_scale,&
25 cp_fm_uplo_to_full
26 USE cp_fm_diag, ONLY: cp_fm_syevd
27 USE cp_fm_types, ONLY: cp_fm_create,&
28 cp_fm_get_info,&
29 cp_fm_release,&
30 cp_fm_set_all,&
31 cp_fm_to_fm,&
32 cp_fm_type
33 USE cp_log_handling, ONLY: cp_get_default_logger,&
34 cp_logger_get_default_unit_nr,&
35 cp_logger_type
36 USE input_constants, ONLY: do_arnoldi,&
37 do_bch,&
38 do_cn,&
39 do_em,&
40 do_etrs,&
41 do_exact,&
42 do_pade,&
43 do_taylor
44 USE iterate_matrix, ONLY: matrix_sqrt_newton_schulz
45 USE kinds, ONLY: dp
46 USE matrix_exp, ONLY: get_nsquare_norder
47 USE parallel_gemm_api, ONLY: parallel_gemm
48 USE qs_environment_types, ONLY: get_qs_env,&
49 qs_environment_type
50 USE qs_mo_types, ONLY: mo_set_type
51 USE rt_make_propagators, ONLY: compute_exponential,&
52 compute_exponential_sparse,&
53 propagate_arnoldi
54 USE rt_propagation_methods, ONLY: calc_sinvh,&
55 put_data_to_history,&
56 s_matrices_create
57 USE rt_propagation_types, ONLY: get_rtp,&
58 rt_prop_release_mos,&
59 rt_prop_type
60#include "../base/base_uses.f90"
61
62 IMPLICIT NONE
63
64 PRIVATE
65
66 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagator_init'
67
68 PUBLIC :: init_propagators, &
69 rt_initialize_rho_from_mos
70
71CONTAINS
72
73! **************************************************************************************************
74!> \brief prepares the initial matrices for the propagators
75!> \param qs_env ...
76!> \author Florian Schiffmann (02.09)
77! **************************************************************************************************
78
79 SUBROUTINE init_propagators(qs_env)
80
81 TYPE(qs_environment_type), POINTER :: qs_env
82
83 INTEGER :: i, imat, unit_nr
84 REAL(kind=dp) :: dt, prefac
85 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new, mos_next, mos_old
86 TYPE(cp_logger_type), POINTER :: logger
87 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: exp_h_new, exp_h_old, matrix_ks, &
88 matrix_ks_im, propagator_matrix, &
89 rho_old, s_mat
90 TYPE(dft_control_type), POINTER :: dft_control
91 TYPE(rt_prop_type), POINTER :: rtp
92 TYPE(rtp_control_type), POINTER :: rtp_control
93
94 CALL get_qs_env(qs_env, &
95 rtp=rtp, &
96 dft_control=dft_control, &
97 matrix_s=s_mat, &
98 matrix_ks=matrix_ks, &
99 matrix_ks_im=matrix_ks_im)
100
101 rtp_control => dft_control%rtp_control
102 CALL get_rtp(rtp=rtp, exp_h_old=exp_h_old, exp_h_new=exp_h_new, &
103 propagator_matrix=propagator_matrix, dt=dt)
104 CALL s_matrices_create(s_mat, rtp)
105 CALL calc_sinvh(rtp, matrix_ks, matrix_ks_im, rtp_control)
106 DO i = 1, SIZE(exp_h_old)
107 CALL dbcsr_copy(exp_h_old(i)%matrix, exp_h_new(i)%matrix)
108 END DO
109 ! use the fact that CN propagator is a first order pade approximation on the EM propagator
110 IF (rtp_control%propagator == do_cn) THEN
111 rtp%orders(1, :) = 0; rtp%orders(2, :) = 1; rtp_control%mat_exp = do_pade; rtp_control%propagator = do_em
112 ELSE IF (rtp_control%mat_exp == do_pade .OR. rtp_control%mat_exp == do_taylor) THEN
113 IF (rtp%linear_scaling) THEN
114 CALL get_maxabs_eigval_sparse(rtp, s_mat, matrix_ks, rtp_control)
115 ELSE
116 CALL get_maxabs_eigval(rtp, s_mat, matrix_ks, rtp_control)
117 END IF
118 END IF
119 ! Safety check for the matrix exponantial scheme wrt the MO representation
120 IF ((rtp_control%mat_exp == do_pade .OR. rtp_control%mat_exp == do_arnoldi) .AND. rtp%linear_scaling) THEN
121 ! get a useful output_unit
122 logger => cp_get_default_logger()
123 IF (logger%para_env%is_source()) THEN
124 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
125 WRITE (unit_nr, *) "linear_scaling currently does not support pade nor arnoldi exponentials, switching to taylor"
126 END IF
127 rtp_control%mat_exp = do_taylor
128
129 ELSE IF (rtp_control%mat_exp == do_bch .AND. .NOT. rtp%linear_scaling) THEN
130 ! get a useful output_unit
131 logger => cp_get_default_logger()
132 IF (logger%para_env%is_source()) THEN
133 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
134 WRITE (unit_nr, *) "BCH exponential currently only support linear_scaling, switching to arnoldi"
135 END IF
136 rtp_control%mat_exp = do_arnoldi
137
138 END IF
139 ! We have no clue yet about next H so we use initial H for t and t+dt
140 ! Due to different nature of the propagator the prefactor has to be adopted
141 SELECT CASE (rtp_control%propagator)
142 CASE (do_etrs)
143 prefac = -0.5_dp*dt
144 CASE (do_em)
145 prefac = -1.0_dp*dt
146 END SELECT
147 DO imat = 1, SIZE(exp_h_new)
148 CALL dbcsr_copy(propagator_matrix(imat)%matrix, exp_h_new(imat)%matrix)
149 CALL dbcsr_scale(propagator_matrix(imat)%matrix, prefac)
150 END DO
151
152 ! For ETRS this bit could be avoided but it drastically simplifies the workflow afterwards.
153 ! If we compute the half propagated mos/exponential already here, we ensure everything is done
154 ! with the correct S matrix and all information as during RTP/EMD are computed.
155 ! Therefore we might accept to compute an unnesscesary exponential but understand the code afterwards
156 IF (rtp_control%propagator == do_etrs) THEN
157 IF (rtp_control%mat_exp == do_arnoldi) THEN
158 rtp%iter = 0
159 CALL propagate_arnoldi(rtp, rtp_control)
160 CALL get_rtp(rtp=rtp, mos_new=mos_new, mos_next=mos_next)
161 DO imat = 1, SIZE(mos_new)
162 CALL cp_fm_to_fm(mos_new(imat), mos_next(imat))
163 END DO
164 ELSEIF (rtp_control%mat_exp == do_bch .OR. rtp_control%mat_exp == do_exact) THEN
165 ELSE
166 IF (rtp%linear_scaling) THEN
167 CALL compute_exponential_sparse(exp_h_new, propagator_matrix, rtp_control, rtp)
168 ELSE
169 CALL compute_exponential(exp_h_new, propagator_matrix, rtp_control, rtp)
170 END IF
171 DO imat = 1, SIZE(exp_h_new)
172 CALL dbcsr_copy(exp_h_old(imat)%matrix, exp_h_new(imat)%matrix)
173 END DO
174 END IF
175 END IF
176
177 IF (rtp%linear_scaling) THEN
178 CALL get_rtp(rtp=rtp, rho_old=rho_old)
179 ELSE
180 CALL get_rtp(rtp=rtp, mos_old=mos_old)
181 END IF
182 CALL put_data_to_history(rtp, mos=mos_old, s_mat=s_mat, ihist=1, rho=rho_old)
183
184 END SUBROUTINE init_propagators
185
186! **************************************************************************************************
187!> \brief gets an estimate for the 2-norm of KS (diagnaliztion of KS) and
188!> calculates the order and number of squaring steps for Taylor or
189!> Pade matrix exponential
190!> \param rtp ...
191!> \param s_mat ...
192!> \param matrix_ks ...
193!> \param rtp_control ...
194!> \author Florian Schiffmann (02.09)
195! **************************************************************************************************
196
197 SUBROUTINE get_maxabs_eigval(rtp, s_mat, matrix_ks, rtp_control)
198 TYPE(rt_prop_type), POINTER :: rtp
199 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: s_mat, matrix_ks
200 TYPE(rtp_control_type), POINTER :: rtp_control
201
202 CHARACTER(len=*), PARAMETER :: routinen = 'get_maxabs_eigval'
203 REAL(kind=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
204
205 INTEGER :: handle, ispin, method, ndim
206 LOGICAL :: emd
207 REAL(dp) :: max_eval, min_eval, norm2, scale, t
208 REAL(dp), ALLOCATABLE, DIMENSION(:) :: eigval_h
209 TYPE(cp_fm_type) :: eigvec_h, h_fm, s_half, s_inv_fm, &
210 s_minus_half, tmp, tmp_mat_h
211 TYPE(dbcsr_type), POINTER :: s_inv
212
213 CALL timeset(routinen, handle)
214
215 CALL get_rtp(rtp=rtp, s_inv=s_inv, dt=t)
216
217 CALL cp_fm_create(s_inv_fm, &
218 matrix_struct=rtp%ao_ao_fmstruct, &
219 name="S_inv")
220 CALL copy_dbcsr_to_fm(s_inv, s_inv_fm)
221
222 CALL cp_fm_create(s_half, &
223 matrix_struct=rtp%ao_ao_fmstruct, &
224 name="S_half")
225
226 CALL cp_fm_create(s_minus_half, &
227 matrix_struct=rtp%ao_ao_fmstruct, &
228 name="S_minus_half")
229
230 CALL cp_fm_create(h_fm, &
231 matrix_struct=rtp%ao_ao_fmstruct, &
232 name="RTP_H_FM")
233
234 CALL cp_fm_create(tmp_mat_h, &
235 matrix_struct=rtp%ao_ao_fmstruct, &
236 name="TMP_H")
237
238 ndim = s_inv_fm%matrix_struct%nrow_global
239 scale = 1.0_dp
240 IF (rtp_control%propagator == do_etrs) scale = 2.0_dp
241 t = -t/scale
242
243 ! Create the overlap matrices
244
245 CALL cp_fm_create(tmp, &
246 matrix_struct=rtp%ao_ao_fmstruct, &
247 name="tmp_mat")
248
249 CALL cp_fm_create(eigvec_h, &
250 matrix_struct=rtp%ao_ao_fmstruct, &
251 name="tmp_EVEC")
252
253 ALLOCATE (eigval_h(ndim))
254 CALL copy_dbcsr_to_fm(s_mat(1)%matrix, tmp)
255 CALL cp_fm_uplo_to_full(tmp, eigvec_h)
256
257 CALL cp_fm_syevd(tmp, eigvec_h, eigval_h)
258
259 eigval_h(:) = one/eigval_h(:)
260 CALL backtransform_matrix(eigval_h, eigvec_h, s_inv_fm)
261 eigval_h(:) = sqrt(eigval_h(:))
262 CALL backtransform_matrix(eigval_h, eigvec_h, s_minus_half)
263 eigval_h(:) = one/eigval_h(:)
264 CALL backtransform_matrix(eigval_h, eigvec_h, s_half)
265 CALL cp_fm_release(eigvec_h)
266 CALL cp_fm_release(tmp)
267
268 IF (rtp_control%mat_exp == do_taylor) method = 1
269 IF (rtp_control%mat_exp == do_pade) method = 2
270 emd = (.NOT. rtp_control%fixed_ions)
271
272 DO ispin = 1, SIZE(matrix_ks)
273
274 CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, h_fm)
275 CALL cp_fm_uplo_to_full(h_fm, tmp_mat_h)
276 CALL cp_fm_scale(t, h_fm)
277
278 CALL parallel_gemm("N", "N", ndim, ndim, ndim, one, h_fm, s_minus_half, zero, &
279 tmp_mat_h)
280 CALL parallel_gemm("N", "N", ndim, ndim, ndim, one, s_minus_half, tmp_mat_h, zero, &
281 h_fm)
282
283 CALL cp_fm_syevd(h_fm, tmp_mat_h, eigval_h)
284 min_eval = minval(eigval_h)
285 max_eval = maxval(eigval_h)
286 norm2 = 2.0_dp*max(abs(min_eval), abs(max_eval))
287 CALL get_nsquare_norder(norm2, rtp%orders(1, ispin), rtp%orders(2, ispin), &
288 rtp_control%eps_exp, method, emd)
289 END DO
290
291 DEALLOCATE (eigval_h)
292
293 CALL copy_fm_to_dbcsr(s_inv_fm, s_inv)
294 CALL cp_fm_release(s_inv_fm)
295 CALL cp_fm_release(s_half)
296 CALL cp_fm_release(s_minus_half)
297 CALL cp_fm_release(h_fm)
298 CALL cp_fm_release(tmp_mat_h)
299
300 CALL timestop(handle)
301
302 END SUBROUTINE get_maxabs_eigval
303
304! **************************************************************************************************
305!> \brief gets an estimate for the 2-norm of KS (diagnaliztion of KS) and
306!> calculates the order and number of squaring steps for Taylor or
307!> Pade matrix exponential. Based on the full matrix code.
308!> \param rtp ...
309!> \param s_mat ...
310!> \param matrix_ks ...
311!> \param rtp_control ...
312!> \author Samuel Andermatt (02.14)
313! **************************************************************************************************
314
315 SUBROUTINE get_maxabs_eigval_sparse(rtp, s_mat, matrix_ks, rtp_control)
316 TYPE(rt_prop_type), POINTER :: rtp
317 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: s_mat, matrix_ks
318 TYPE(rtp_control_type), POINTER :: rtp_control
319
320 CHARACTER(len=*), PARAMETER :: routinen = 'get_maxabs_eigval_sparse'
321 REAL(kind=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
322
323 INTEGER :: handle, ispin, method
324 LOGICAL :: converged, emd
325 REAL(dp) :: max_ev, min_ev, norm2, scale, t
326 TYPE(dbcsr_type), POINTER :: s_half, s_minus_half, tmp, tmp2
327
328 CALL timeset(routinen, handle)
329
330 CALL get_rtp(rtp=rtp, dt=t)
331
332 NULLIFY (s_half)
333 ALLOCATE (s_half)
334 CALL dbcsr_create(s_half, template=s_mat(1)%matrix)
335 NULLIFY (s_minus_half)
336 ALLOCATE (s_minus_half)
337 CALL dbcsr_create(s_minus_half, template=s_mat(1)%matrix)
338 NULLIFY (tmp)
339 ALLOCATE (tmp)
340 CALL dbcsr_create(tmp, template=s_mat(1)%matrix, matrix_type="N")
341 NULLIFY (tmp2)
342 ALLOCATE (tmp2)
343 CALL dbcsr_create(tmp2, template=s_mat(1)%matrix, matrix_type="N")
344 scale = 1.0_dp
345 IF (rtp_control%propagator == do_etrs) scale = 2.0_dp
346 t = -t/scale
347 emd = (.NOT. rtp_control%fixed_ions)
348
349 IF (rtp_control%mat_exp == do_taylor) method = 1
350 IF (rtp_control%mat_exp == do_pade) method = 2
351 CALL matrix_sqrt_newton_schulz(s_half, s_minus_half, s_mat(1)%matrix, rtp%filter_eps, &
352 rtp%newton_schulz_order, rtp%lanzcos_threshold, rtp%lanzcos_max_iter)
353 DO ispin = 1, SIZE(matrix_ks)
354 CALL dbcsr_multiply("N", "N", t, matrix_ks(ispin)%matrix, s_minus_half, zero, tmp, &
355 filter_eps=rtp%filter_eps)
356 CALL dbcsr_multiply("N", "N", one, s_minus_half, tmp, zero, tmp2, &
357 filter_eps=rtp%filter_eps)
358 CALL arnoldi_extremal(tmp2, max_ev, min_ev, threshold=rtp%lanzcos_threshold, &
359 max_iter=rtp%lanzcos_max_iter, converged=converged)
360 norm2 = 2.0_dp*max(abs(min_ev), abs(max_ev))
361 CALL get_nsquare_norder(norm2, rtp%orders(1, ispin), rtp%orders(2, ispin), &
362 rtp_control%eps_exp, method, emd)
363 END DO
364
365 CALL dbcsr_deallocate_matrix(s_half)
366 CALL dbcsr_deallocate_matrix(s_minus_half)
367 CALL dbcsr_deallocate_matrix(tmp)
368 CALL dbcsr_deallocate_matrix(tmp2)
369
370 CALL timestop(handle)
371
372 END SUBROUTINE get_maxabs_eigval_sparse
373
374! **************************************************************************************************
375!> \brief Is still left from diagonalization, should be removed later but is
376!> still needed for the guess for the pade/Taylor method
377!> \param Eval ...
378!> \param eigenvec ...
379!> \param matrix ...
380!> \author Florian Schiffmann (02.09)
381! **************************************************************************************************
382
383 SUBROUTINE backtransform_matrix(Eval, eigenvec, matrix)
384
385 REAL(dp), DIMENSION(:), INTENT(in) :: eval
386 TYPE(cp_fm_type), INTENT(IN) :: eigenvec, matrix
387
388 CHARACTER(len=*), PARAMETER :: routinen = 'backtransform_matrix'
389 REAL(kind=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
390
391 INTEGER :: handle, i, j, l, ncol_local, ndim, &
392 nrow_local
393 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
394 TYPE(cp_fm_type) :: tmp
395
396 CALL timeset(routinen, handle)
397 CALL cp_fm_create(tmp, &
398 matrix_struct=matrix%matrix_struct, &
399 name="TMP_BT")
400 CALL cp_fm_get_info(matrix, nrow_local=nrow_local, ncol_local=ncol_local, &
401 row_indices=row_indices, col_indices=col_indices)
402
403 ndim = matrix%matrix_struct%nrow_global
404
405 CALL cp_fm_set_all(tmp, zero, zero)
406 DO i = 1, ncol_local
407 l = col_indices(i)
408 DO j = 1, nrow_local
409 tmp%local_data(j, i) = eigenvec%local_data(j, i)*eval(l)
410 END DO
411 END DO
412 CALL parallel_gemm("N", "T", ndim, ndim, ndim, one, tmp, eigenvec, zero, &
413 matrix)
414
415 CALL cp_fm_release(tmp)
416 CALL timestop(handle)
417
418 END SUBROUTINE backtransform_matrix
419
420! **************************************************************************************************
421!> \brief Computes the density matrix from the mos
422!> Update: Initialized the density matrix from complex mos (for
423!> instance after delta kick)
424!> \param rtp ...
425!> \param mos ...
426!> \param mos_old ...
427!> \author Samuel Andermatt (08.15)
428!> \author Guillaume Le Breton (01.23)
429! **************************************************************************************************
430
431 SUBROUTINE rt_initialize_rho_from_mos(rtp, mos, mos_old)
432
433 TYPE(rt_prop_type), POINTER :: rtp
434 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
435 TYPE(cp_fm_type), DIMENSION(:), OPTIONAL, POINTER :: mos_old
436
437 INTEGER :: im, ispin, ncol, re
438 REAL(kind=dp) :: alpha
439 TYPE(cp_fm_type) :: mos_occ, mos_occ_im
440 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new, rho_old
441
442 CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
443
444 IF (PRESENT(mos_old)) THEN
445 ! Used the mos from delta kick. Initialize both real and im part
446 DO ispin = 1, SIZE(mos_old)/2
447 re = 2*ispin - 1; im = 2*ispin
448 CALL dbcsr_set(rho_old(re)%matrix, 0.0_dp)
449 CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=ncol)
450 CALL cp_fm_create(mos_occ, &
451 matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
452 name="mos_occ")
453 !Real part of rho
454 CALL cp_fm_to_fm(mos_old(re), mos_occ)
455 IF (mos(ispin)%uniform_occupation) THEN
456 alpha = 3.0_dp - real(SIZE(mos_old)/2, dp)
457 CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
458 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
459 matrix_v=mos_occ, &
460 ncol=ncol, &
461 alpha=alpha, keep_sparsity=.false.)
462 ELSE
463 alpha = 1.0_dp
464 CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
465 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
466 matrix_v=mos_old(re), &
467 matrix_g=mos_occ, &
468 ncol=ncol, &
469 alpha=alpha, keep_sparsity=.false.)
470 END IF
471
472 ! Add complex part of MOs, i*i=-1
473 CALL cp_fm_to_fm(mos_old(im), mos_occ)
474 IF (mos(ispin)%uniform_occupation) THEN
475 alpha = 3.0_dp - real(SIZE(mos_old)/2, dp)
476 CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
477 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
478 matrix_v=mos_occ, &
479 ncol=ncol, &
480 alpha=alpha, keep_sparsity=.false.)
481 ELSE
482 alpha = 1.0_dp
483 CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
484 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
485 matrix_v=mos_old(im), &
486 matrix_g=mos_occ, &
487 ncol=ncol, &
488 alpha=alpha, keep_sparsity=.false.)
489 END IF
490 CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
491 CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
492
493 ! Imaginary part of rho
494 CALL cp_fm_create(mos_occ_im, &
495 matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
496 name="mos_occ_im")
497 alpha = 3.0_dp - real(SIZE(mos_old)/2, dp)
498 CALL cp_fm_to_fm(mos_old(re), mos_occ)
499 CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
500 CALL cp_fm_to_fm(mos_old(im), mos_occ_im)
501 CALL cp_fm_column_scale(mos_occ_im, mos(ispin)%occupation_numbers/alpha)
502 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(im)%matrix, &
503 matrix_v=mos_occ_im, &
504 matrix_g=mos_occ, &
505 ncol=ncol, &
506 alpha=2.0_dp*alpha, &
507 symmetry_mode=-1, keep_sparsity=.false.)
508
509 CALL dbcsr_filter(rho_old(im)%matrix, rtp%filter_eps)
510 CALL dbcsr_copy(rho_new(im)%matrix, rho_old(im)%matrix)
511 CALL cp_fm_release(mos_occ_im)
512 CALL cp_fm_release(mos_occ)
513 END DO
514 ! Release the mos used to apply the delta kick, no longer required
515 CALL rt_prop_release_mos(rtp)
516 ELSE
517 DO ispin = 1, SIZE(mos)
518 re = 2*ispin - 1
519 CALL dbcsr_set(rho_old(re)%matrix, 0.0_dp)
520 CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=ncol)
521
522 CALL cp_fm_create(mos_occ, &
523 matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
524 name="mos_occ")
525 CALL cp_fm_to_fm(mos(ispin)%mo_coeff, mos_occ)
526 IF (mos(ispin)%uniform_occupation) THEN
527 alpha = 3.0_dp - real(SIZE(mos), dp)
528 CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
529 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
530 matrix_v=mos_occ, &
531 ncol=ncol, &
532 alpha=alpha, keep_sparsity=.false.)
533 ELSE
534 alpha = 1.0_dp
535 CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
536 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
537 matrix_v=mos(ispin)%mo_coeff, &
538 matrix_g=mos_occ, &
539 ncol=ncol, &
540 alpha=alpha, keep_sparsity=.false.)
541 END IF
542 CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
543 CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
544 CALL cp_fm_release(mos_occ)
545 END DO
546 END IF
547
548 END SUBROUTINE rt_initialize_rho_from_mos
549
550END MODULE rt_propagator_init
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
arnoldi_api
arnoldi iteration using dbcsr
Definition arnoldi_api.F:16
arnoldi_api::arnoldi_extremal
subroutine, public arnoldi_extremal(matrix_a, max_ev, min_ev, converged, threshold, max_iter)
simple wrapper to estimate extremal eigenvalues with arnoldi, using the old lanczos interface this hi...
Definition arnoldi_api.F:254
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1178
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:370
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1086
cp_dbcsr_api::dbcsr_filter
subroutine, public dbcsr_filter(matrix, eps)
...
Definition cp_dbcsr_api.F:657
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1194
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:214
cp_dbcsr_operations::cp_dbcsr_plus_fm_fm_t
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
Definition cp_dbcsr_operations.F:695
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:111
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition cp_fm_basic_linalg.F:1937
cp_fm_basic_linalg::cp_fm_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1801
cp_fm_basic_linalg::cp_fm_scale
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
Definition cp_fm_basic_linalg.F:1719
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::cp_fm_syevd
subroutine, public cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
Computes all eigenvalues and vectors of a real symmetric matrix significantly faster than syevx,...
Definition cp_fm_diag.F:443
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:166
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1064
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:533
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:567
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_bch
integer, parameter, public do_bch
Definition input_constants.F:928
input_constants::do_etrs
integer, parameter, public do_etrs
Definition input_constants.F:934
input_constants::do_pade
integer, parameter, public do_pade
Definition input_constants.F:928
input_constants::do_cn
integer, parameter, public do_cn
Definition input_constants.F:934
input_constants::do_exact
integer, parameter, public do_exact
Definition input_constants.F:928
input_constants::do_taylor
integer, parameter, public do_taylor
Definition input_constants.F:928
input_constants::do_em
integer, parameter, public do_em
Definition input_constants.F:934
input_constants::do_arnoldi
integer, parameter, public do_arnoldi
Definition input_constants.F:928
iterate_matrix
Routines useful for iterative matrix calculations.
Definition iterate_matrix.F:13
iterate_matrix::matrix_sqrt_newton_schulz
subroutine, public matrix_sqrt_newton_schulz(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, eps_lanczos, max_iter_lanczos, symmetrize, converged, iounit)
compute the sqrt of a matrix via the sign function and the corresponding Newton-Schulz iterations the...
Definition iterate_matrix.F:1626
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
matrix_exp
Routines for calculating a complex matrix exponential.
Definition matrix_exp.F:13
matrix_exp::get_nsquare_norder
subroutine, public get_nsquare_norder(norm, nsquare, norder, eps_exp, method, do_emd)
optimization function for pade/taylor order and number of squaring steps
Definition matrix_exp.F:244
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:526
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
rt_make_propagators
Routines for calculating a complex matrix exponential.
Definition rt_make_propagators.F:13
rt_make_propagators::compute_exponential_sparse
subroutine, public compute_exponential_sparse(propagator, propagator_matrix, rtp_control, rtp)
Sparse versions of the matrix exponentials.
Definition rt_make_propagators.F:403
rt_make_propagators::compute_exponential
subroutine, public compute_exponential(propagator, propagator_matrix, rtp_control, rtp)
decides which type of exponential has to be computed
Definition rt_make_propagators.F:339
rt_make_propagators::propagate_arnoldi
subroutine, public propagate_arnoldi(rtp, rtp_control)
computes U_prop*MOs using arnoldi subspace algorithm
Definition rt_make_propagators.F:200
rt_propagation_methods
Routines for propagating the orbitals.
Definition rt_propagation_methods.F:12
rt_propagation_methods::calc_sinvh
subroutine, public calc_sinvh(rtp, matrix_ks, matrix_ks_im, rtp_control)
computes S_inv*H, if needed Sinv*B and S_inv*H_imag and store these quantities to the
Definition rt_propagation_methods.F:407
rt_propagation_methods::put_data_to_history
subroutine, public put_data_to_history(rtp, mos, rho, s_mat, ihist)
...
Definition rt_propagation_methods.F:829
rt_propagation_methods::s_matrices_create
subroutine, public s_matrices_create(s_mat, rtp)
calculates the needed overlap-like matrices depending on the way the exponential is calculated,...
Definition rt_propagation_methods.F:489
rt_propagation_types
Types and set_get for real time propagation depending on runtype and diagonalization method different...
Definition rt_propagation_types.F:21
rt_propagation_types::rt_prop_release_mos
subroutine, public rt_prop_release_mos(rtp)
Deallocated the mos for rtp...
Definition rt_propagation_types.F:515
rt_propagation_types::get_rtp
subroutine, public get_rtp(rtp, exp_h_old, exp_h_new, h_last_iter, rho_old, rho_next, rho_new, mos, mos_new, mos_old, mos_next, s_inv, s_half, s_minus_half, b_mat, c_mat, propagator_matrix, mixing, mixing_factor, s_der, dt, nsteps, sinvh, sinvh_imag, sinvb, admm_mos)
...
Definition rt_propagation_types.F:416
rt_propagator_init
Routines for that prepare rtp and EMD.
Definition rt_propagator_init.F:12
rt_propagator_init::init_propagators
subroutine, public init_propagators(qs_env)
prepares the initial matrices for the propagators
Definition rt_propagator_init.F:80
rt_propagator_init::rt_initialize_rho_from_mos
subroutine, public rt_initialize_rho_from_mos(rtp, mos, mos_old)
Computes the density matrix from the mos Update: Initialized the density matrix from complex mos (for...
Definition rt_propagator_init.F:432
cp_control_types::dft_control_type
Definition cp_control_types.F:624
cp_control_types::rtp_control_type
Definition cp_control_types.F:91
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47
rt_propagation_types::rt_prop_type
Definition rt_propagation_types.F:76