(git:374b731)
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-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! **************************************************************************************************
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_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
70CONTAINS
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
549END MODULE rt_propagator_init
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
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:1882
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:1759
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:1680
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
cp_control_types::dft_control_type
Definition cp_control_types.F:559
cp_control_types::rtp_control_type
Definition cp_control_types.F:65
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
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:215
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
rt_propagation_types::rt_prop_type
Definition rt_propagation_types.F:76