(git:e7e05ae)
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-2024 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 ! **************************************************************************************************
12 MODULE 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_operations, ONLY: copy_dbcsr_to_fm,&
18  copy_fm_to_dbcsr,&
19  cp_dbcsr_plus_fm_fm_t
20  USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
21  cp_fm_scale,&
22  cp_fm_upper_to_full
23  USE cp_fm_diag, ONLY: cp_fm_syevd
24  USE cp_fm_types, ONLY: cp_fm_create,&
25  cp_fm_get_info,&
26  cp_fm_release,&
27  cp_fm_set_all,&
28  cp_fm_to_fm,&
29  cp_fm_type
30  USE cp_log_handling, ONLY: cp_get_default_logger,&
31  cp_logger_get_default_unit_nr,&
32  cp_logger_type
33  USE dbcsr_api, ONLY: &
34  dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_filter, dbcsr_multiply, &
35  dbcsr_p_type, dbcsr_scale, dbcsr_set, dbcsr_type
36  USE input_constants, ONLY: do_arnoldi,&
37  do_bch,&
38  do_cn,&
39  do_em,&
40  do_etrs,&
41  do_pade,&
42  do_taylor
43  USE iterate_matrix, ONLY: matrix_sqrt_newton_schulz
44  USE kinds, ONLY: dp
45  USE matrix_exp, ONLY: get_nsquare_norder
46  USE parallel_gemm_api, ONLY: parallel_gemm
47  USE qs_environment_types, ONLY: get_qs_env,&
48  qs_environment_type
49  USE qs_mo_types, ONLY: mo_set_type
50  USE rt_make_propagators, ONLY: compute_exponential,&
51  compute_exponential_sparse,&
52  propagate_arnoldi
53  USE rt_propagation_methods, ONLY: calc_sinvh,&
54  put_data_to_history,&
55  s_matrices_create
56  USE rt_propagation_types, ONLY: get_rtp,&
57  rt_prop_release_mos,&
58  rt_prop_type
59 #include "../base/base_uses.f90"
60 
61  IMPLICIT NONE
62 
63  PRIVATE
64 
65  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagator_init'
66 
67  PUBLIC :: init_propagators, &
68  rt_initialize_rho_from_mos
69 
70 CONTAINS
71 
72 ! **************************************************************************************************
73 !> \brief prepares the initial matrices for the propagators
74 !> \param qs_env ...
75 !> \author Florian Schiffmann (02.09)
76 ! **************************************************************************************************
77 
78  SUBROUTINE init_propagators(qs_env)
79 
80  TYPE(qs_environment_type), POINTER :: qs_env
81 
82  INTEGER :: i, imat, unit_nr
83  REAL(kind=dp) :: dt, prefac
84  TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new, mos_next, mos_old
85  TYPE(cp_logger_type), POINTER :: logger
86  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: exp_h_new, exp_h_old, matrix_ks, &
87  matrix_ks_im, propagator_matrix, &
88  rho_old, s_mat
89  TYPE(dft_control_type), POINTER :: dft_control
90  TYPE(rt_prop_type), POINTER :: rtp
91  TYPE(rtp_control_type), POINTER :: rtp_control
92 
93  CALL get_qs_env(qs_env, &
94  rtp=rtp, &
95  dft_control=dft_control, &
96  matrix_s=s_mat, &
97  matrix_ks=matrix_ks, &
98  matrix_ks_im=matrix_ks_im)
99 
100  rtp_control => dft_control%rtp_control
101  CALL get_rtp(rtp=rtp, exp_h_old=exp_h_old, exp_h_new=exp_h_new, &
102  propagator_matrix=propagator_matrix, dt=dt)
103  CALL s_matrices_create(s_mat, rtp)
104  CALL calc_sinvh(rtp, matrix_ks, matrix_ks_im, rtp_control)
105  DO i = 1, SIZE(exp_h_old)
106  CALL dbcsr_copy(exp_h_old(i)%matrix, exp_h_new(i)%matrix)
107  END DO
108  ! use the fact that CN propagator is a first order pade approximation on the EM propagator
109  IF (rtp_control%propagator == do_cn) THEN
110  rtp%orders(1, :) = 0; rtp%orders(2, :) = 1; rtp_control%mat_exp = do_pade; rtp_control%propagator = do_em
111  ELSE IF (rtp_control%mat_exp == do_pade .OR. rtp_control%mat_exp == do_taylor) THEN
112  IF (rtp%linear_scaling) THEN
113  CALL get_maxabs_eigval_sparse(rtp, s_mat, matrix_ks, rtp_control)
114  ELSE
115  CALL get_maxabs_eigval(rtp, s_mat, matrix_ks, rtp_control)
116  END IF
117  END IF
118  ! Safety check for the matrix exponantial scheme wrt the MO representation
119  IF ((rtp_control%mat_exp == do_pade .OR. rtp_control%mat_exp == do_arnoldi) .AND. rtp%linear_scaling) THEN
120  ! get a useful output_unit
121  logger => cp_get_default_logger()
122  IF (logger%para_env%is_source()) THEN
123  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
124  WRITE (unit_nr, *) "linear_scaling currently does not support pade nor arnoldi exponentials, switching to taylor"
125  END IF
126  rtp_control%mat_exp = do_taylor
127 
128  ELSE IF (rtp_control%mat_exp == do_bch .AND. .NOT. rtp%linear_scaling) THEN
129  ! get a useful output_unit
130  logger => cp_get_default_logger()
131  IF (logger%para_env%is_source()) THEN
132  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
133  WRITE (unit_nr, *) "BCH exponential currently only support linear_scaling, switching to arnoldi"
134  END IF
135  rtp_control%mat_exp = do_arnoldi
136 
137  END IF
138  ! We have no clue yet about next H so we use initial H for t and t+dt
139  ! Due to different nature of the propagator the prefactor has to be adopted
140  SELECT CASE (rtp_control%propagator)
141  CASE (do_etrs)
142  prefac = -0.5_dp*dt
143  CASE (do_em)
144  prefac = -1.0_dp*dt
145  END SELECT
146  DO imat = 1, SIZE(exp_h_new)
147  CALL dbcsr_copy(propagator_matrix(imat)%matrix, exp_h_new(imat)%matrix)
148  CALL dbcsr_scale(propagator_matrix(imat)%matrix, prefac)
149  END DO
150 
151  ! For ETRS this bit could be avoided but it drastically simplifies the workflow afterwards.
152  ! If we compute the half propagated mos/exponential already here, we ensure everything is done
153  ! with the correct S matrix and all information as during RTP/EMD are computed.
154  ! Therefore we might accept to compute an unnesscesary exponential but understand the code afterwards
155  IF (rtp_control%propagator == do_etrs) THEN
156  IF (rtp_control%mat_exp == do_arnoldi) THEN
157  rtp%iter = 0
158  CALL propagate_arnoldi(rtp, rtp_control)
159  CALL get_rtp(rtp=rtp, mos_new=mos_new, mos_next=mos_next)
160  DO imat = 1, SIZE(mos_new)
161  CALL cp_fm_to_fm(mos_new(imat), mos_next(imat))
162  END DO
163  ELSEIF (rtp_control%mat_exp == do_bch) THEN
164  ELSE
165  IF (rtp%linear_scaling) THEN
166  CALL compute_exponential_sparse(exp_h_new, propagator_matrix, rtp_control, rtp)
167  ELSE
168  CALL compute_exponential(exp_h_new, propagator_matrix, rtp_control, rtp)
169  END IF
170  DO imat = 1, SIZE(exp_h_new)
171  CALL dbcsr_copy(exp_h_old(imat)%matrix, exp_h_new(imat)%matrix)
172  END DO
173  END IF
174  END IF
175 
176  IF (rtp%linear_scaling) THEN
177  CALL get_rtp(rtp=rtp, rho_old=rho_old)
178  ELSE
179  CALL get_rtp(rtp=rtp, mos_old=mos_old)
180  END IF
181  CALL put_data_to_history(rtp, mos=mos_old, s_mat=s_mat, ihist=1, rho=rho_old)
182 
183  END SUBROUTINE init_propagators
184 
185 ! **************************************************************************************************
186 !> \brief gets an estimate for the 2-norm of KS (diagnaliztion of KS) and
187 !> calculates the order and number of squaring steps for Taylor or
188 !> Pade matrix exponential
189 !> \param rtp ...
190 !> \param s_mat ...
191 !> \param matrix_ks ...
192 !> \param rtp_control ...
193 !> \author Florian Schiffmann (02.09)
194 ! **************************************************************************************************
195 
196  SUBROUTINE get_maxabs_eigval(rtp, s_mat, matrix_ks, rtp_control)
197  TYPE(rt_prop_type), POINTER :: rtp
198  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: s_mat, matrix_ks
199  TYPE(rtp_control_type), POINTER :: rtp_control
200 
201  CHARACTER(len=*), PARAMETER :: routinen = 'get_maxabs_eigval'
202  REAL(kind=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
203 
204  INTEGER :: handle, ispin, method, ndim
205  LOGICAL :: emd
206  REAL(dp) :: max_eval, min_eval, norm2, scale, t
207  REAL(dp), ALLOCATABLE, DIMENSION(:) :: eigval_h
208  TYPE(cp_fm_type) :: eigvec_h, h_fm, s_half, s_inv_fm, &
209  s_minus_half, tmp, tmp_mat_h
210  TYPE(dbcsr_type), POINTER :: s_inv
211 
212  CALL timeset(routinen, handle)
213 
214  CALL get_rtp(rtp=rtp, s_inv=s_inv, dt=t)
215 
216  CALL cp_fm_create(s_inv_fm, &
217  matrix_struct=rtp%ao_ao_fmstruct, &
218  name="S_inv")
219  CALL copy_dbcsr_to_fm(s_inv, s_inv_fm)
220 
221  CALL cp_fm_create(s_half, &
222  matrix_struct=rtp%ao_ao_fmstruct, &
223  name="S_half")
224 
225  CALL cp_fm_create(s_minus_half, &
226  matrix_struct=rtp%ao_ao_fmstruct, &
227  name="S_minus_half")
228 
229  CALL cp_fm_create(h_fm, &
230  matrix_struct=rtp%ao_ao_fmstruct, &
231  name="RTP_H_FM")
232 
233  CALL cp_fm_create(tmp_mat_h, &
234  matrix_struct=rtp%ao_ao_fmstruct, &
235  name="TMP_H")
236 
237  ndim = s_inv_fm%matrix_struct%nrow_global
238  scale = 1.0_dp
239  IF (rtp_control%propagator == do_etrs) scale = 2.0_dp
240  t = -t/scale
241 
242  ! Create the overlap matrices
243 
244  CALL cp_fm_create(tmp, &
245  matrix_struct=rtp%ao_ao_fmstruct, &
246  name="tmp_mat")
247 
248  CALL cp_fm_create(eigvec_h, &
249  matrix_struct=rtp%ao_ao_fmstruct, &
250  name="tmp_EVEC")
251 
252  ALLOCATE (eigval_h(ndim))
253  CALL copy_dbcsr_to_fm(s_mat(1)%matrix, tmp)
254  CALL cp_fm_upper_to_full(tmp, eigvec_h)
255 
256  CALL cp_fm_syevd(tmp, eigvec_h, eigval_h)
257 
258  eigval_h(:) = one/eigval_h(:)
259  CALL backtransform_matrix(eigval_h, eigvec_h, s_inv_fm)
260  eigval_h(:) = sqrt(eigval_h(:))
261  CALL backtransform_matrix(eigval_h, eigvec_h, s_minus_half)
262  eigval_h(:) = one/eigval_h(:)
263  CALL backtransform_matrix(eigval_h, eigvec_h, s_half)
264  CALL cp_fm_release(eigvec_h)
265  CALL cp_fm_release(tmp)
266 
267  IF (rtp_control%mat_exp == do_taylor) method = 1
268  IF (rtp_control%mat_exp == do_pade) method = 2
269  emd = (.NOT. rtp_control%fixed_ions)
270 
271  DO ispin = 1, SIZE(matrix_ks)
272 
273  CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, h_fm)
274  CALL cp_fm_upper_to_full(h_fm, tmp_mat_h)
275  CALL cp_fm_scale(t, h_fm)
276 
277  CALL parallel_gemm("N", "N", ndim, ndim, ndim, one, h_fm, s_minus_half, zero, &
278  tmp_mat_h)
279  CALL parallel_gemm("N", "N", ndim, ndim, ndim, one, s_minus_half, tmp_mat_h, zero, &
280  h_fm)
281 
282  CALL cp_fm_syevd(h_fm, tmp_mat_h, eigval_h)
283  min_eval = minval(eigval_h)
284  max_eval = maxval(eigval_h)
285  norm2 = 2.0_dp*max(abs(min_eval), abs(max_eval))
286  CALL get_nsquare_norder(norm2, rtp%orders(1, ispin), rtp%orders(2, ispin), &
287  rtp_control%eps_exp, method, emd)
288  END DO
289 
290  DEALLOCATE (eigval_h)
291 
292  CALL copy_fm_to_dbcsr(s_inv_fm, s_inv)
293  CALL cp_fm_release(s_inv_fm)
294  CALL cp_fm_release(s_half)
295  CALL cp_fm_release(s_minus_half)
296  CALL cp_fm_release(h_fm)
297  CALL cp_fm_release(tmp_mat_h)
298 
299  CALL timestop(handle)
300 
301  END SUBROUTINE get_maxabs_eigval
302 
303 ! **************************************************************************************************
304 !> \brief gets an estimate for the 2-norm of KS (diagnaliztion of KS) and
305 !> calculates the order and number of squaring steps for Taylor or
306 !> Pade matrix exponential. Based on the full matrix code.
307 !> \param rtp ...
308 !> \param s_mat ...
309 !> \param matrix_ks ...
310 !> \param rtp_control ...
311 !> \author Samuel Andermatt (02.14)
312 ! **************************************************************************************************
313 
314  SUBROUTINE get_maxabs_eigval_sparse(rtp, s_mat, matrix_ks, rtp_control)
315  TYPE(rt_prop_type), POINTER :: rtp
316  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: s_mat, matrix_ks
317  TYPE(rtp_control_type), POINTER :: rtp_control
318 
319  CHARACTER(len=*), PARAMETER :: routinen = 'get_maxabs_eigval_sparse'
320  REAL(kind=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
321 
322  INTEGER :: handle, ispin, method
323  LOGICAL :: converged, emd
324  REAL(dp) :: max_ev, min_ev, norm2, scale, t
325  TYPE(dbcsr_type), POINTER :: s_half, s_minus_half, tmp, tmp2
326 
327  CALL timeset(routinen, handle)
328 
329  CALL get_rtp(rtp=rtp, dt=t)
330 
331  NULLIFY (s_half)
332  ALLOCATE (s_half)
333  CALL dbcsr_create(s_half, template=s_mat(1)%matrix)
334  NULLIFY (s_minus_half)
335  ALLOCATE (s_minus_half)
336  CALL dbcsr_create(s_minus_half, template=s_mat(1)%matrix)
337  NULLIFY (tmp)
338  ALLOCATE (tmp)
339  CALL dbcsr_create(tmp, template=s_mat(1)%matrix, matrix_type="N")
340  NULLIFY (tmp2)
341  ALLOCATE (tmp2)
342  CALL dbcsr_create(tmp2, template=s_mat(1)%matrix, matrix_type="N")
343  scale = 1.0_dp
344  IF (rtp_control%propagator == do_etrs) scale = 2.0_dp
345  t = -t/scale
346  emd = (.NOT. rtp_control%fixed_ions)
347 
348  IF (rtp_control%mat_exp == do_taylor) method = 1
349  IF (rtp_control%mat_exp == do_pade) method = 2
350  CALL matrix_sqrt_newton_schulz(s_half, s_minus_half, s_mat(1)%matrix, rtp%filter_eps, &
351  rtp%newton_schulz_order, rtp%lanzcos_threshold, rtp%lanzcos_max_iter)
352  DO ispin = 1, SIZE(matrix_ks)
353  CALL dbcsr_multiply("N", "N", t, matrix_ks(ispin)%matrix, s_minus_half, zero, tmp, &
354  filter_eps=rtp%filter_eps)
355  CALL dbcsr_multiply("N", "N", one, s_minus_half, tmp, zero, tmp2, &
356  filter_eps=rtp%filter_eps)
357  CALL arnoldi_extremal(tmp2, max_ev, min_ev, threshold=rtp%lanzcos_threshold, &
358  max_iter=rtp%lanzcos_max_iter, converged=converged)
359  norm2 = 2.0_dp*max(abs(min_ev), abs(max_ev))
360  CALL get_nsquare_norder(norm2, rtp%orders(1, ispin), rtp%orders(2, ispin), &
361  rtp_control%eps_exp, method, emd)
362  END DO
363 
364  CALL dbcsr_deallocate_matrix(s_half)
365  CALL dbcsr_deallocate_matrix(s_minus_half)
366  CALL dbcsr_deallocate_matrix(tmp)
367  CALL dbcsr_deallocate_matrix(tmp2)
368 
369  CALL timestop(handle)
370 
371  END SUBROUTINE
372 
373 ! **************************************************************************************************
374 !> \brief Is still left from diagonalization, should be removed later but is
375 !> still needed for the guess for the pade/Taylor method
376 !> \param Eval ...
377 !> \param eigenvec ...
378 !> \param matrix ...
379 !> \author Florian Schiffmann (02.09)
380 ! **************************************************************************************************
381 
382  SUBROUTINE backtransform_matrix(Eval, eigenvec, matrix)
383 
384  REAL(dp), DIMENSION(:), INTENT(in) :: eval
385  TYPE(cp_fm_type), INTENT(IN) :: eigenvec, matrix
386 
387  CHARACTER(len=*), PARAMETER :: routinen = 'backtransform_matrix'
388  REAL(kind=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
389 
390  INTEGER :: handle, i, j, l, ncol_local, ndim, &
391  nrow_local
392  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
393  TYPE(cp_fm_type) :: tmp
394 
395  CALL timeset(routinen, handle)
396  CALL cp_fm_create(tmp, &
397  matrix_struct=matrix%matrix_struct, &
398  name="TMP_BT")
399  CALL cp_fm_get_info(matrix, nrow_local=nrow_local, ncol_local=ncol_local, &
400  row_indices=row_indices, col_indices=col_indices)
401 
402  ndim = matrix%matrix_struct%nrow_global
403 
404  CALL cp_fm_set_all(tmp, zero, zero)
405  DO i = 1, ncol_local
406  l = col_indices(i)
407  DO j = 1, nrow_local
408  tmp%local_data(j, i) = eigenvec%local_data(j, i)*eval(l)
409  END DO
410  END DO
411  CALL parallel_gemm("N", "T", ndim, ndim, ndim, one, tmp, eigenvec, zero, &
412  matrix)
413 
414  CALL cp_fm_release(tmp)
415  CALL timestop(handle)
416 
417  END SUBROUTINE backtransform_matrix
418 
419 ! **************************************************************************************************
420 !> \brief Computes the density matrix from the mos
421 !> Update: Initialized the density matrix from complex mos (for
422 !> instance after delta kick)
423 !> \param rtp ...
424 !> \param mos ...
425 !> \param mos_old ...
426 !> \author Samuel Andermatt (08.15)
427 !> \author Guillaume Le Breton (01.23)
428 ! **************************************************************************************************
429 
430  SUBROUTINE rt_initialize_rho_from_mos(rtp, mos, mos_old)
431 
432  TYPE(rt_prop_type), POINTER :: rtp
433  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
434  TYPE(cp_fm_type), DIMENSION(:), OPTIONAL, POINTER :: mos_old
435 
436  INTEGER :: im, ispin, ncol, re
437  REAL(kind=dp) :: alpha
438  TYPE(cp_fm_type) :: mos_occ, mos_occ_im
439  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new, rho_old
440 
441  CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
442 
443  IF (PRESENT(mos_old)) THEN
444  ! Used the mos from delta kick. Initialize both real and im part
445  DO ispin = 1, SIZE(mos_old)/2
446  re = 2*ispin - 1; im = 2*ispin
447  CALL dbcsr_set(rho_old(re)%matrix, 0.0_dp)
448  CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=ncol)
449  CALL cp_fm_create(mos_occ, &
450  matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
451  name="mos_occ")
452  !Real part of rho
453  CALL cp_fm_to_fm(mos_old(re), mos_occ)
454  IF (mos(ispin)%uniform_occupation) THEN
455  alpha = 3.0_dp - real(SIZE(mos_old)/2, dp)
456  CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
457  CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
458  matrix_v=mos_occ, &
459  ncol=ncol, &
460  alpha=alpha, keep_sparsity=.false.)
461  ELSE
462  alpha = 1.0_dp
463  CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
464  CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
465  matrix_v=mos_old(re), &
466  matrix_g=mos_occ, &
467  ncol=ncol, &
468  alpha=alpha, keep_sparsity=.false.)
469  END IF
470 
471  ! Add complex part of MOs, i*i=-1
472  CALL cp_fm_to_fm(mos_old(im), mos_occ)
473  IF (mos(ispin)%uniform_occupation) THEN
474  alpha = 3.0_dp - real(SIZE(mos_old)/2, dp)
475  CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
476  CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
477  matrix_v=mos_occ, &
478  ncol=ncol, &
479  alpha=alpha, keep_sparsity=.false.)
480  ELSE
481  alpha = 1.0_dp
482  CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
483  CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
484  matrix_v=mos_old(im), &
485  matrix_g=mos_occ, &
486  ncol=ncol, &
487  alpha=alpha, keep_sparsity=.false.)
488  END IF
489  CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
490  CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
491 
492  ! Imaginary part of rho
493  CALL cp_fm_create(mos_occ_im, &
494  matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
495  name="mos_occ_im")
496  alpha = 3.0_dp - real(SIZE(mos_old)/2, dp)
497  CALL cp_fm_to_fm(mos_old(re), mos_occ)
498  CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
499  CALL cp_fm_to_fm(mos_old(im), mos_occ_im)
500  CALL cp_fm_column_scale(mos_occ_im, mos(ispin)%occupation_numbers/alpha)
501  CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(im)%matrix, &
502  matrix_v=mos_occ_im, &
503  matrix_g=mos_occ, &
504  ncol=ncol, &
505  alpha=2.0_dp*alpha, &
506  symmetry_mode=-1, keep_sparsity=.false.)
507 
508  CALL dbcsr_filter(rho_old(im)%matrix, rtp%filter_eps)
509  CALL dbcsr_copy(rho_new(im)%matrix, rho_old(im)%matrix)
510  CALL cp_fm_release(mos_occ_im)
511  CALL cp_fm_release(mos_occ)
512  END DO
513  ! Release the mos used to apply the delta kick, no longer required
514  CALL rt_prop_release_mos(rtp)
515  ELSE
516  DO ispin = 1, SIZE(mos)
517  re = 2*ispin - 1
518  CALL dbcsr_set(rho_old(re)%matrix, 0.0_dp)
519  CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=ncol)
520 
521  CALL cp_fm_create(mos_occ, &
522  matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
523  name="mos_occ")
524  CALL cp_fm_to_fm(mos(ispin)%mo_coeff, mos_occ)
525  IF (mos(ispin)%uniform_occupation) THEN
526  alpha = 3.0_dp - real(SIZE(mos), dp)
527  CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
528  CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
529  matrix_v=mos_occ, &
530  ncol=ncol, &
531  alpha=alpha, keep_sparsity=.false.)
532  ELSE
533  alpha = 1.0_dp
534  CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
535  CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
536  matrix_v=mos(ispin)%mo_coeff, &
537  matrix_g=mos_occ, &
538  ncol=ncol, &
539  alpha=alpha, keep_sparsity=.false.)
540  END IF
541  CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
542  CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
543  CALL cp_fm_release(mos_occ)
544  END DO
545  END IF
546 
547  END SUBROUTINE rt_initialize_rho_from_mos
548 
549 END MODULE rt_propagator_init
arnoldi_api
arnoldi iteration using dbcsr
Definition: arnoldi_api.F:15
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_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:208
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:887
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:118
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:1816
cp_fm_basic_linalg::cp_fm_upper_to_full
subroutine, public cp_fm_upper_to_full(matrix, work)
given an upper triangular matrix computes the corresponding full matrix
Definition: cp_fm_basic_linalg.F:1693
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:1614
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:413
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
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:1016
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:535
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
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:918
input_constants::do_etrs
integer, parameter, public do_etrs
Definition: input_constants.F:923
input_constants::do_pade
integer, parameter, public do_pade
Definition: input_constants.F:918
input_constants::do_cn
integer, parameter, public do_cn
Definition: input_constants.F:923
input_constants::do_taylor
integer, parameter, public do_taylor
Definition: input_constants.F:918
input_constants::do_em
integer, parameter, public do_em
Definition: input_constants.F:923
input_constants::do_arnoldi
integer, parameter, public do_arnoldi
Definition: input_constants.F:918
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)
compute the sqrt of a matrix via the sign function and the corresponding Newton-Schulz iterations the...
Definition: iterate_matrix.F:1624
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_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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, ecoul_1c, rho0_s_rs, rho0_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, 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, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
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:406
rt_propagation_methods::put_data_to_history
subroutine, public put_data_to_history(rtp, mos, rho, s_mat, ihist)
...
Definition: rt_propagation_methods.F:828
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:488
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:507
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:408
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:79
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:431