Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines for 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_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 :
71 : CONTAINS
72 :
73 : ! **************************************************************************************************
74 : !> \brief prepares the initial matrices for the propagators
75 : !> \param qs_env ...
76 : !> \author Florian Schiffmann (02.09)
77 : ! **************************************************************************************************
78 :
79 396 : 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 198 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new, mos_next, mos_old
86 : TYPE(cp_logger_type), POINTER :: logger
87 198 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: exp_H_new, exp_H_old, matrix_ks, &
88 198 : matrix_ks_im, propagator_matrix, &
89 198 : 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 198 : matrix_ks_im=matrix_ks_im)
100 :
101 198 : 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 198 : propagator_matrix=propagator_matrix, dt=dt)
104 198 : CALL s_matrices_create(s_mat, rtp)
105 198 : CALL calc_SinvH(rtp, matrix_ks, matrix_ks_im, rtp_control)
106 730 : DO i = 1, SIZE(exp_H_old)
107 730 : 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 198 : IF (rtp_control%propagator == do_cn) THEN
111 12 : rtp%orders(1, :) = 0; rtp%orders(2, :) = 1; rtp_control%mat_exp = do_pade; rtp_control%propagator = do_em
112 194 : ELSE IF (rtp_control%mat_exp == do_pade .OR. rtp_control%mat_exp == do_taylor) THEN
113 154 : IF (rtp%linear_scaling) THEN
114 76 : CALL get_maxabs_eigval_sparse(rtp, s_mat, matrix_ks, rtp_control)
115 : ELSE
116 78 : 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 198 : 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 0 : logger => cp_get_default_logger()
123 0 : IF (logger%para_env%is_source()) THEN
124 0 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
125 0 : WRITE (unit_nr, *) "linear_scaling currently does not support pade nor arnoldi exponentials, switching to taylor"
126 : END IF
127 0 : rtp_control%mat_exp = do_taylor
128 :
129 198 : ELSE IF (rtp_control%mat_exp == do_bch .AND. .NOT. rtp%linear_scaling) THEN
130 : ! get a useful output_unit
131 0 : logger => cp_get_default_logger()
132 0 : IF (logger%para_env%is_source()) THEN
133 0 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
134 0 : WRITE (unit_nr, *) "BCH exponential currently only support linear_scaling, switching to arnoldi"
135 : END IF
136 0 : 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 390 : SELECT CASE (rtp_control%propagator)
142 : CASE (do_etrs)
143 192 : prefac = -0.5_dp*dt
144 : CASE (do_em)
145 198 : prefac = -1.0_dp*dt
146 : END SELECT
147 730 : DO imat = 1, SIZE(exp_H_new)
148 532 : CALL dbcsr_copy(propagator_matrix(imat)%matrix, exp_H_new(imat)%matrix)
149 730 : 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 198 : IF (rtp_control%propagator == do_etrs) THEN
157 192 : IF (rtp_control%mat_exp == do_arnoldi) THEN
158 26 : rtp%iter = 0
159 26 : CALL propagate_arnoldi(rtp, rtp_control)
160 26 : CALL get_rtp(rtp=rtp, mos_new=mos_new, mos_next=mos_next)
161 98 : DO imat = 1, SIZE(mos_new)
162 98 : CALL cp_fm_to_fm(mos_new(imat), mos_next(imat))
163 : END DO
164 166 : ELSEIF (rtp_control%mat_exp == do_bch .OR. rtp_control%mat_exp == do_exact) THEN
165 : ELSE
166 152 : IF (rtp%linear_scaling) THEN
167 76 : CALL compute_exponential_sparse(exp_H_new, propagator_matrix, rtp_control, rtp)
168 : ELSE
169 76 : CALL compute_exponential(exp_H_new, propagator_matrix, rtp_control, rtp)
170 : END IF
171 568 : DO imat = 1, SIZE(exp_H_new)
172 568 : CALL dbcsr_copy(exp_H_old(imat)%matrix, exp_H_new(imat)%matrix)
173 : END DO
174 : END IF
175 : END IF
176 :
177 198 : IF (rtp%linear_scaling) THEN
178 90 : CALL get_rtp(rtp=rtp, rho_old=rho_old)
179 : ELSE
180 108 : CALL get_rtp(rtp=rtp, mos_old=mos_old)
181 : END IF
182 198 : CALL put_data_to_history(rtp, mos=mos_old, s_mat=s_mat, ihist=1, rho=rho_old)
183 :
184 198 : 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 78 : 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 78 : CALL timeset(routineN, handle)
214 :
215 78 : 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 78 : name="S_inv")
220 78 : 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 78 : name="S_half")
225 :
226 : CALL cp_fm_create(S_minus_half, &
227 : matrix_struct=rtp%ao_ao_fmstruct, &
228 78 : name="S_minus_half")
229 :
230 : CALL cp_fm_create(H_fm, &
231 : matrix_struct=rtp%ao_ao_fmstruct, &
232 78 : name="RTP_H_FM")
233 :
234 : CALL cp_fm_create(tmp_mat_H, &
235 : matrix_struct=rtp%ao_ao_fmstruct, &
236 78 : name="TMP_H")
237 :
238 78 : ndim = S_inv_fm%matrix_struct%nrow_global
239 78 : scale = 1.0_dp
240 78 : IF (rtp_control%propagator == do_etrs) scale = 2.0_dp
241 78 : t = -t/scale
242 :
243 : ! Create the overlap matrices
244 :
245 : CALL cp_fm_create(tmp, &
246 : matrix_struct=rtp%ao_ao_fmstruct, &
247 78 : name="tmp_mat")
248 :
249 : CALL cp_fm_create(eigvec_H, &
250 : matrix_struct=rtp%ao_ao_fmstruct, &
251 78 : name="tmp_EVEC")
252 :
253 234 : ALLOCATE (eigval_H(ndim))
254 78 : CALL copy_dbcsr_to_fm(s_mat(1)%matrix, tmp)
255 78 : CALL cp_fm_uplo_to_full(tmp, eigvec_H)
256 :
257 78 : CALL cp_fm_syevd(tmp, eigvec_H, eigval_H)
258 :
259 1150 : eigval_H(:) = one/eigval_H(:)
260 78 : CALL backtransform_matrix(eigval_H, eigvec_H, S_inv_fm)
261 1150 : eigval_H(:) = SQRT(eigval_H(:))
262 78 : CALL backtransform_matrix(eigval_H, eigvec_H, S_minus_half)
263 1150 : eigval_H(:) = one/eigval_H(:)
264 78 : CALL backtransform_matrix(eigval_H, eigvec_H, S_half)
265 78 : CALL cp_fm_release(eigvec_H)
266 78 : CALL cp_fm_release(tmp)
267 :
268 78 : IF (rtp_control%mat_exp == do_taylor) method = 1
269 78 : IF (rtp_control%mat_exp == do_pade) method = 2
270 78 : emd = (.NOT. rtp_control%fixed_ions)
271 :
272 176 : DO ispin = 1, SIZE(matrix_ks)
273 :
274 98 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, H_fm)
275 98 : CALL cp_fm_uplo_to_full(H_fm, tmp_mat_H)
276 98 : 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 98 : tmp_mat_H)
280 : CALL parallel_gemm("N", "N", ndim, ndim, ndim, one, S_minus_half, tmp_mat_H, zero, &
281 98 : H_fm)
282 :
283 98 : CALL cp_fm_syevd(H_fm, tmp_mat_H, eigval_H)
284 1628 : min_eval = MINVAL(eigval_H)
285 1628 : max_eval = MAXVAL(eigval_H)
286 98 : 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 176 : rtp_control%eps_exp, method, emd)
289 : END DO
290 :
291 78 : DEALLOCATE (eigval_H)
292 :
293 78 : CALL copy_fm_to_dbcsr(S_inv_fm, S_inv)
294 78 : CALL cp_fm_release(S_inv_fm)
295 78 : CALL cp_fm_release(S_half)
296 78 : CALL cp_fm_release(S_minus_half)
297 78 : CALL cp_fm_release(H_fm)
298 78 : CALL cp_fm_release(tmp_mat_H)
299 :
300 78 : CALL timestop(handle)
301 :
302 234 : 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 152 : 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 76 : CALL timeset(routineN, handle)
329 :
330 76 : CALL get_rtp(rtp=rtp, dt=t)
331 :
332 : NULLIFY (s_half)
333 76 : ALLOCATE (s_half)
334 76 : CALL dbcsr_create(s_half, template=s_mat(1)%matrix)
335 : NULLIFY (s_minus_half)
336 76 : ALLOCATE (s_minus_half)
337 76 : CALL dbcsr_create(s_minus_half, template=s_mat(1)%matrix)
338 : NULLIFY (tmp)
339 76 : ALLOCATE (tmp)
340 76 : CALL dbcsr_create(tmp, template=s_mat(1)%matrix, matrix_type="N")
341 : NULLIFY (tmp2)
342 76 : ALLOCATE (tmp2)
343 76 : CALL dbcsr_create(tmp2, template=s_mat(1)%matrix, matrix_type="N")
344 76 : scale = 1.0_dp
345 76 : IF (rtp_control%propagator == do_etrs) scale = 2.0_dp
346 76 : t = -t/scale
347 76 : emd = (.NOT. rtp_control%fixed_ions)
348 :
349 76 : IF (rtp_control%mat_exp == do_taylor) method = 1
350 76 : 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 76 : rtp%newton_schulz_order, rtp%lanzcos_threshold, rtp%lanzcos_max_iter)
353 188 : DO ispin = 1, SIZE(matrix_ks)
354 : CALL dbcsr_multiply("N", "N", t, matrix_ks(ispin)%matrix, s_minus_half, zero, tmp, &
355 112 : filter_eps=rtp%filter_eps)
356 : CALL dbcsr_multiply("N", "N", one, s_minus_half, tmp, zero, tmp2, &
357 112 : filter_eps=rtp%filter_eps)
358 : CALL arnoldi_extremal(tmp2, max_ev, min_ev, threshold=rtp%lanzcos_threshold, &
359 112 : max_iter=rtp%lanzcos_max_iter, converged=converged)
360 112 : 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 188 : rtp_control%eps_exp, method, emd)
363 : END DO
364 :
365 76 : CALL dbcsr_deallocate_matrix(s_half)
366 76 : CALL dbcsr_deallocate_matrix(s_minus_half)
367 76 : CALL dbcsr_deallocate_matrix(tmp)
368 76 : CALL dbcsr_deallocate_matrix(tmp2)
369 :
370 76 : CALL timestop(handle)
371 :
372 76 : END SUBROUTINE
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 234 : 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 234 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
394 : TYPE(cp_fm_type) :: tmp
395 :
396 234 : CALL timeset(routineN, handle)
397 : CALL cp_fm_create(tmp, &
398 : matrix_struct=matrix%matrix_struct, &
399 234 : name="TMP_BT")
400 : CALL cp_fm_get_info(matrix, nrow_local=nrow_local, ncol_local=ncol_local, &
401 234 : row_indices=row_indices, col_indices=col_indices)
402 :
403 234 : ndim = matrix%matrix_struct%nrow_global
404 :
405 234 : CALL cp_fm_set_all(tmp, zero, zero)
406 3450 : DO i = 1, ncol_local
407 3216 : l = col_indices(i)
408 30276 : DO j = 1, nrow_local
409 30042 : 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 234 : matrix)
414 :
415 234 : CALL cp_fm_release(tmp)
416 234 : CALL timestop(handle)
417 :
418 234 : 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 64 : 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 64 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new, rho_old
441 :
442 64 : CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
443 :
444 64 : IF (PRESENT(mos_old)) THEN
445 : ! Used the mos from delta kick. Initialize both real and im part
446 106 : DO ispin = 1, SIZE(mos_old)/2
447 66 : re = 2*ispin - 1; im = 2*ispin
448 66 : CALL dbcsr_set(rho_old(re)%matrix, 0.0_dp)
449 66 : 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 66 : name="mos_occ")
453 : !Real part of rho
454 66 : CALL cp_fm_to_fm(mos_old(re), mos_occ)
455 66 : IF (mos(ispin)%uniform_occupation) THEN
456 58 : alpha = 3.0_dp - REAL(SIZE(mos_old)/2, dp)
457 280 : 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 58 : alpha=alpha, keep_sparsity=.FALSE.)
462 : ELSE
463 8 : alpha = 1.0_dp
464 116 : 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 8 : alpha=alpha, keep_sparsity=.FALSE.)
470 : END IF
471 :
472 : ! Add complex part of MOs, i*i=-1
473 66 : CALL cp_fm_to_fm(mos_old(im), mos_occ)
474 66 : IF (mos(ispin)%uniform_occupation) THEN
475 58 : alpha = 3.0_dp - REAL(SIZE(mos_old)/2, dp)
476 280 : 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 58 : alpha=alpha, keep_sparsity=.FALSE.)
481 : ELSE
482 8 : alpha = 1.0_dp
483 116 : 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 8 : alpha=alpha, keep_sparsity=.FALSE.)
489 : END IF
490 66 : CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
491 66 : 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 66 : name="mos_occ_im")
497 66 : alpha = 3.0_dp - REAL(SIZE(mos_old)/2, dp)
498 66 : CALL cp_fm_to_fm(mos_old(re), mos_occ)
499 396 : CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
500 66 : CALL cp_fm_to_fm(mos_old(im), mos_occ_im)
501 396 : 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 66 : symmetry_mode=-1, keep_sparsity=.FALSE.)
508 :
509 66 : CALL dbcsr_filter(rho_old(im)%matrix, rtp%filter_eps)
510 66 : CALL dbcsr_copy(rho_new(im)%matrix, rho_old(im)%matrix)
511 66 : CALL cp_fm_release(mos_occ_im)
512 238 : CALL cp_fm_release(mos_occ)
513 : END DO
514 : ! Release the mos used to apply the delta kick, no longer required
515 40 : CALL rt_prop_release_mos(rtp)
516 : ELSE
517 50 : DO ispin = 1, SIZE(mos)
518 26 : re = 2*ispin - 1
519 26 : CALL dbcsr_set(rho_old(re)%matrix, 0.0_dp)
520 26 : 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 26 : name="mos_occ")
525 26 : CALL cp_fm_to_fm(mos(ispin)%mo_coeff, mos_occ)
526 26 : IF (mos(ispin)%uniform_occupation) THEN
527 22 : alpha = 3.0_dp - REAL(SIZE(mos), dp)
528 110 : 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 22 : alpha=alpha, keep_sparsity=.FALSE.)
533 : ELSE
534 4 : alpha = 1.0_dp
535 88 : 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 4 : alpha=alpha, keep_sparsity=.FALSE.)
541 : END IF
542 26 : CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
543 26 : CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
544 76 : CALL cp_fm_release(mos_occ)
545 : END DO
546 : END IF
547 :
548 64 : END SUBROUTINE rt_initialize_rho_from_mos
549 :
550 : END MODULE rt_propagator_init
|
