Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2023 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 356 : 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 178 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new, mos_next, mos_old
85 : TYPE(cp_logger_type), POINTER :: logger
86 178 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: exp_H_new, exp_H_old, matrix_ks, &
87 178 : matrix_ks_im, propagator_matrix, &
88 178 : 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 178 : matrix_ks_im=matrix_ks_im)
99 :
100 178 : 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 178 : propagator_matrix=propagator_matrix, dt=dt)
103 178 : CALL s_matrices_create(s_mat, rtp)
104 178 : CALL calc_SinvH(rtp, matrix_ks, matrix_ks_im, rtp_control)
105 654 : DO i = 1, SIZE(exp_H_old)
106 654 : 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 178 : IF (rtp_control%propagator == do_cn) THEN
110 12 : rtp%orders(1, :) = 0; rtp%orders(2, :) = 1; rtp_control%mat_exp = do_pade; rtp_control%propagator = do_em
111 174 : ELSE IF (rtp_control%mat_exp == do_pade .OR. rtp_control%mat_exp == do_taylor) THEN
112 142 : IF (rtp%linear_scaling) THEN
113 76 : CALL get_maxabs_eigval_sparse(rtp, s_mat, matrix_ks, rtp_control)
114 : ELSE
115 66 : CALL get_maxabs_eigval(rtp, s_mat, matrix_ks, rtp_control)
116 : END IF
117 : END IF
118 178 : IF (rtp_control%mat_exp == do_pade .AND. rtp%linear_scaling) THEN
119 : ! get a useful output_unit
120 0 : logger => cp_get_default_logger()
121 0 : IF (logger%para_env%is_source()) THEN
122 0 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
123 0 : WRITE (unit_nr, *) "linear_scaling currently does not support pade exponentials, switching to taylor"
124 : END IF
125 0 : rtp_control%mat_exp = do_taylor
126 : END IF
127 : ! We have no clue yet about next H so we use initial H for t and t+dt
128 : ! Due to different nature of the propagator the prefactor has to be adopted
129 350 : SELECT CASE (rtp_control%propagator)
130 : CASE (do_etrs)
131 172 : prefac = -0.5_dp*dt
132 : CASE (do_em)
133 178 : prefac = -1.0_dp*dt
134 : END SELECT
135 654 : DO imat = 1, SIZE(exp_H_new)
136 476 : CALL dbcsr_copy(propagator_matrix(imat)%matrix, exp_H_new(imat)%matrix)
137 654 : CALL dbcsr_scale(propagator_matrix(imat)%matrix, prefac)
138 : END DO
139 :
140 : ! For ETRS this bit could be avoided but it drastically simplifies the workflow afterwards.
141 : ! If we compute the half propagated mos/exponential already here, we ensure everything is done
142 : ! with the correct S matrix and all information as during RTP/EMD are computed.
143 : ! Therefore we might accept to compute an unnesscesary exponential but understand the code afterwards
144 178 : IF (rtp_control%propagator == do_etrs) THEN
145 172 : IF (rtp_control%mat_exp == do_arnoldi) THEN
146 18 : rtp%iter = 0
147 18 : CALL propagate_arnoldi(rtp, rtp_control)
148 18 : CALL get_rtp(rtp=rtp, mos_new=mos_new, mos_next=mos_next)
149 66 : DO imat = 1, SIZE(mos_new)
150 66 : CALL cp_fm_to_fm(mos_new(imat), mos_next(imat))
151 : END DO
152 154 : ELSEIF (rtp_control%mat_exp == do_bch) THEN
153 : ELSE
154 140 : IF (rtp%linear_scaling) THEN
155 76 : CALL compute_exponential_sparse(exp_H_new, propagator_matrix, rtp_control, rtp)
156 : ELSE
157 64 : CALL compute_exponential(exp_H_new, propagator_matrix, rtp_control, rtp)
158 : END IF
159 524 : DO imat = 1, SIZE(exp_H_new)
160 524 : CALL dbcsr_copy(exp_H_old(imat)%matrix, exp_H_new(imat)%matrix)
161 : END DO
162 : END IF
163 : END IF
164 :
165 178 : IF (rtp%linear_scaling) THEN
166 90 : CALL get_rtp(rtp=rtp, rho_old=rho_old)
167 : ELSE
168 88 : CALL get_rtp(rtp=rtp, mos_old=mos_old)
169 : END IF
170 178 : CALL put_data_to_history(rtp, mos=mos_old, s_mat=s_mat, ihist=1, rho=rho_old)
171 :
172 178 : END SUBROUTINE init_propagators
173 :
174 : ! **************************************************************************************************
175 : !> \brief gets an estimate for the 2-norm of KS (diagnaliztion of KS) and
176 : !> calculates the order and number of squaring steps for Taylor or
177 : !> Pade matrix exponential
178 : !> \param rtp ...
179 : !> \param s_mat ...
180 : !> \param matrix_ks ...
181 : !> \param rtp_control ...
182 : !> \author Florian Schiffmann (02.09)
183 : ! **************************************************************************************************
184 :
185 66 : SUBROUTINE get_maxabs_eigval(rtp, s_mat, matrix_ks, rtp_control)
186 : TYPE(rt_prop_type), POINTER :: rtp
187 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: s_mat, matrix_ks
188 : TYPE(rtp_control_type), POINTER :: rtp_control
189 :
190 : CHARACTER(len=*), PARAMETER :: routineN = 'get_maxabs_eigval'
191 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
192 :
193 : INTEGER :: handle, ispin, method, ndim
194 : LOGICAL :: emd
195 : REAL(dp) :: max_eval, min_eval, norm2, scale, t
196 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: eigval_H
197 : TYPE(cp_fm_type) :: eigvec_H, H_fm, S_half, S_inv_fm, &
198 : S_minus_half, tmp, tmp_mat_H
199 : TYPE(dbcsr_type), POINTER :: S_inv
200 :
201 66 : CALL timeset(routineN, handle)
202 :
203 66 : CALL get_rtp(rtp=rtp, S_inv=S_inv, dt=t)
204 :
205 : CALL cp_fm_create(S_inv_fm, &
206 : matrix_struct=rtp%ao_ao_fmstruct, &
207 66 : name="S_inv")
208 66 : CALL copy_dbcsr_to_fm(S_inv, S_inv_fm)
209 :
210 : CALL cp_fm_create(S_half, &
211 : matrix_struct=rtp%ao_ao_fmstruct, &
212 66 : name="S_half")
213 :
214 : CALL cp_fm_create(S_minus_half, &
215 : matrix_struct=rtp%ao_ao_fmstruct, &
216 66 : name="S_minus_half")
217 :
218 : CALL cp_fm_create(H_fm, &
219 : matrix_struct=rtp%ao_ao_fmstruct, &
220 66 : name="RTP_H_FM")
221 :
222 : CALL cp_fm_create(tmp_mat_H, &
223 : matrix_struct=rtp%ao_ao_fmstruct, &
224 66 : name="TMP_H")
225 :
226 66 : ndim = S_inv_fm%matrix_struct%nrow_global
227 66 : scale = 1.0_dp
228 66 : IF (rtp_control%propagator == do_etrs) scale = 2.0_dp
229 66 : t = -t/scale
230 :
231 : ! Create the overlap matrices
232 :
233 : CALL cp_fm_create(tmp, &
234 : matrix_struct=rtp%ao_ao_fmstruct, &
235 66 : name="tmp_mat")
236 :
237 : CALL cp_fm_create(eigvec_H, &
238 : matrix_struct=rtp%ao_ao_fmstruct, &
239 66 : name="tmp_EVEC")
240 :
241 198 : ALLOCATE (eigval_H(ndim))
242 66 : CALL copy_dbcsr_to_fm(s_mat(1)%matrix, tmp)
243 66 : CALL cp_fm_upper_to_full(tmp, eigvec_H)
244 :
245 66 : CALL cp_fm_syevd(tmp, eigvec_H, eigval_H)
246 :
247 962 : eigval_H(:) = one/eigval_H(:)
248 66 : CALL backtransform_matrix(eigval_H, eigvec_H, S_inv_fm)
249 962 : eigval_H(:) = SQRT(eigval_H(:))
250 66 : CALL backtransform_matrix(eigval_H, eigvec_H, S_minus_half)
251 962 : eigval_H(:) = one/eigval_H(:)
252 66 : CALL backtransform_matrix(eigval_H, eigvec_H, S_half)
253 66 : CALL cp_fm_release(eigvec_H)
254 66 : CALL cp_fm_release(tmp)
255 :
256 66 : IF (rtp_control%mat_exp == do_taylor) method = 1
257 66 : IF (rtp_control%mat_exp == do_pade) method = 2
258 66 : emd = (.NOT. rtp_control%fixed_ions)
259 :
260 148 : DO ispin = 1, SIZE(matrix_ks)
261 :
262 82 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, H_fm)
263 82 : CALL cp_fm_upper_to_full(H_fm, tmp_mat_H)
264 82 : CALL cp_fm_scale(t, H_fm)
265 :
266 : CALL parallel_gemm("N", "N", ndim, ndim, ndim, one, H_fm, S_minus_half, zero, &
267 82 : tmp_mat_H)
268 : CALL parallel_gemm("N", "N", ndim, ndim, ndim, one, S_minus_half, tmp_mat_H, zero, &
269 82 : H_fm)
270 :
271 82 : CALL cp_fm_syevd(H_fm, tmp_mat_H, eigval_H)
272 1324 : min_eval = MINVAL(eigval_H)
273 1324 : max_eval = MAXVAL(eigval_H)
274 82 : norm2 = 2.0_dp*MAX(ABS(min_eval), ABS(max_eval))
275 : CALL get_nsquare_norder(norm2, rtp%orders(1, ispin), rtp%orders(2, ispin), &
276 148 : rtp_control%eps_exp, method, emd)
277 : END DO
278 :
279 66 : DEALLOCATE (eigval_H)
280 :
281 66 : CALL copy_fm_to_dbcsr(S_inv_fm, S_inv)
282 66 : CALL cp_fm_release(S_inv_fm)
283 66 : CALL cp_fm_release(S_half)
284 66 : CALL cp_fm_release(S_minus_half)
285 66 : CALL cp_fm_release(H_fm)
286 66 : CALL cp_fm_release(tmp_mat_H)
287 :
288 66 : CALL timestop(handle)
289 :
290 66 : END SUBROUTINE get_maxabs_eigval
291 :
292 : ! **************************************************************************************************
293 : !> \brief gets an estimate for the 2-norm of KS (diagnaliztion of KS) and
294 : !> calculates the order and number of squaring steps for Taylor or
295 : !> Pade matrix exponential. Based on the full matrix code.
296 : !> \param rtp ...
297 : !> \param s_mat ...
298 : !> \param matrix_ks ...
299 : !> \param rtp_control ...
300 : !> \author Samuel Andermatt (02.14)
301 : ! **************************************************************************************************
302 :
303 152 : SUBROUTINE get_maxabs_eigval_sparse(rtp, s_mat, matrix_ks, rtp_control)
304 : TYPE(rt_prop_type), POINTER :: rtp
305 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: s_mat, matrix_ks
306 : TYPE(rtp_control_type), POINTER :: rtp_control
307 :
308 : CHARACTER(len=*), PARAMETER :: routineN = 'get_maxabs_eigval_sparse'
309 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
310 :
311 : INTEGER :: handle, ispin, method
312 : LOGICAL :: converged, emd
313 : REAL(dp) :: max_ev, min_ev, norm2, scale, t
314 : TYPE(dbcsr_type), POINTER :: s_half, s_minus_half, tmp, tmp2
315 :
316 76 : CALL timeset(routineN, handle)
317 :
318 76 : CALL get_rtp(rtp=rtp, dt=t)
319 :
320 : NULLIFY (s_half)
321 76 : ALLOCATE (s_half)
322 76 : CALL dbcsr_create(s_half, template=s_mat(1)%matrix)
323 : NULLIFY (s_minus_half)
324 76 : ALLOCATE (s_minus_half)
325 76 : CALL dbcsr_create(s_minus_half, template=s_mat(1)%matrix)
326 : NULLIFY (tmp)
327 76 : ALLOCATE (tmp)
328 76 : CALL dbcsr_create(tmp, template=s_mat(1)%matrix, matrix_type="N")
329 : NULLIFY (tmp2)
330 76 : ALLOCATE (tmp2)
331 76 : CALL dbcsr_create(tmp2, template=s_mat(1)%matrix, matrix_type="N")
332 76 : scale = 1.0_dp
333 76 : IF (rtp_control%propagator == do_etrs) scale = 2.0_dp
334 76 : t = -t/scale
335 76 : emd = (.NOT. rtp_control%fixed_ions)
336 :
337 76 : IF (rtp_control%mat_exp == do_taylor) method = 1
338 76 : IF (rtp_control%mat_exp == do_pade) method = 2
339 : CALL matrix_sqrt_Newton_Schulz(s_half, s_minus_half, s_mat(1)%matrix, rtp%filter_eps, &
340 76 : rtp%newton_schulz_order, rtp%lanzcos_threshold, rtp%lanzcos_max_iter)
341 188 : DO ispin = 1, SIZE(matrix_ks)
342 : CALL dbcsr_multiply("N", "N", t, matrix_ks(ispin)%matrix, s_minus_half, zero, tmp, &
343 112 : filter_eps=rtp%filter_eps)
344 : CALL dbcsr_multiply("N", "N", one, s_minus_half, tmp, zero, tmp2, &
345 112 : filter_eps=rtp%filter_eps)
346 : CALL arnoldi_extremal(tmp2, max_ev, min_ev, threshold=rtp%lanzcos_threshold, &
347 112 : max_iter=rtp%lanzcos_max_iter, converged=converged)
348 112 : norm2 = 2.0_dp*MAX(ABS(min_ev), ABS(max_ev))
349 : CALL get_nsquare_norder(norm2, rtp%orders(1, ispin), rtp%orders(2, ispin), &
350 188 : rtp_control%eps_exp, method, emd)
351 : END DO
352 :
353 76 : CALL dbcsr_deallocate_matrix(s_half)
354 76 : CALL dbcsr_deallocate_matrix(s_minus_half)
355 76 : CALL dbcsr_deallocate_matrix(tmp)
356 76 : CALL dbcsr_deallocate_matrix(tmp2)
357 :
358 76 : CALL timestop(handle)
359 :
360 76 : END SUBROUTINE
361 :
362 : ! **************************************************************************************************
363 : !> \brief Is still left from diagonalization, should be removed later but is
364 : !> still needed for the guess for the pade/Taylor method
365 : !> \param Eval ...
366 : !> \param eigenvec ...
367 : !> \param matrix ...
368 : !> \author Florian Schiffmann (02.09)
369 : ! **************************************************************************************************
370 :
371 198 : SUBROUTINE backtransform_matrix(Eval, eigenvec, matrix)
372 :
373 : REAL(dp), DIMENSION(:), INTENT(in) :: Eval
374 : TYPE(cp_fm_type), INTENT(IN) :: eigenvec, matrix
375 :
376 : CHARACTER(len=*), PARAMETER :: routineN = 'backtransform_matrix'
377 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
378 :
379 : INTEGER :: handle, i, j, l, ncol_local, ndim, &
380 : nrow_local
381 198 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
382 : TYPE(cp_fm_type) :: tmp
383 :
384 198 : CALL timeset(routineN, handle)
385 : CALL cp_fm_create(tmp, &
386 : matrix_struct=matrix%matrix_struct, &
387 198 : name="TMP_BT")
388 : CALL cp_fm_get_info(matrix, nrow_local=nrow_local, ncol_local=ncol_local, &
389 198 : row_indices=row_indices, col_indices=col_indices)
390 :
391 198 : ndim = matrix%matrix_struct%nrow_global
392 :
393 198 : CALL cp_fm_set_all(tmp, zero, zero)
394 2886 : DO i = 1, ncol_local
395 2688 : l = col_indices(i)
396 25056 : DO j = 1, nrow_local
397 24858 : tmp%local_data(j, i) = eigenvec%local_data(j, i)*Eval(l)
398 : END DO
399 : END DO
400 : CALL parallel_gemm("N", "T", ndim, ndim, ndim, one, tmp, eigenvec, zero, &
401 198 : matrix)
402 :
403 198 : CALL cp_fm_release(tmp)
404 198 : CALL timestop(handle)
405 :
406 198 : END SUBROUTINE backtransform_matrix
407 :
408 : ! **************************************************************************************************
409 : !> \brief Computes the density matrix from the mos
410 : !> Update: Initialized the density matrix from complex mos (for
411 : !> instance after delta kick)
412 : !> \param rtp ...
413 : !> \param mos ...
414 : !> \param mos_old ...
415 : !> \author Samuel Andermatt (08.15)
416 : !> \author Guillaume Le Breton (01.23)
417 : ! **************************************************************************************************
418 :
419 64 : SUBROUTINE rt_initialize_rho_from_mos(rtp, mos, mos_old)
420 :
421 : TYPE(rt_prop_type), POINTER :: rtp
422 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
423 : TYPE(cp_fm_type), DIMENSION(:), OPTIONAL, POINTER :: mos_old
424 :
425 : INTEGER :: im, ispin, ncol, re
426 : REAL(KIND=dp) :: alpha
427 : TYPE(cp_fm_type) :: mos_occ, mos_occ_im
428 64 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new, rho_old
429 :
430 64 : CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
431 :
432 64 : IF (PRESENT(mos_old)) THEN
433 : ! Used the mos from delta kick. Initialize both real and im part
434 106 : DO ispin = 1, SIZE(mos_old)/2
435 66 : re = 2*ispin - 1; im = 2*ispin
436 66 : CALL dbcsr_set(rho_old(re)%matrix, 0.0_dp)
437 66 : CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=ncol)
438 : CALL cp_fm_create(mos_occ, &
439 : matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
440 66 : name="mos_occ")
441 : !Real part of rho
442 66 : CALL cp_fm_to_fm(mos_old(re), mos_occ)
443 66 : IF (mos(ispin)%uniform_occupation) THEN
444 58 : alpha = 3.0_dp - REAL(SIZE(mos_old)/2, dp)
445 280 : CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
446 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
447 : matrix_v=mos_occ, &
448 : ncol=ncol, &
449 58 : alpha=alpha, keep_sparsity=.FALSE.)
450 : ELSE
451 8 : alpha = 1.0_dp
452 116 : CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
453 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
454 : matrix_v=mos_old(re), &
455 : matrix_g=mos_occ, &
456 : ncol=ncol, &
457 8 : alpha=alpha, keep_sparsity=.FALSE.)
458 : END IF
459 :
460 : ! Add complex part of MOs, i*i=-1
461 66 : CALL cp_fm_to_fm(mos_old(im), mos_occ)
462 66 : IF (mos(ispin)%uniform_occupation) THEN
463 58 : alpha = 3.0_dp - REAL(SIZE(mos_old)/2, dp)
464 280 : 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_occ, &
467 : ncol=ncol, &
468 58 : alpha=alpha, keep_sparsity=.FALSE.)
469 : ELSE
470 8 : alpha = 1.0_dp
471 116 : CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
472 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
473 : matrix_v=mos_old(im), &
474 : matrix_g=mos_occ, &
475 : ncol=ncol, &
476 8 : alpha=alpha, keep_sparsity=.FALSE.)
477 : END IF
478 66 : CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
479 66 : CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
480 :
481 : ! Imaginary part of rho
482 : CALL cp_fm_create(mos_occ_im, &
483 : matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
484 66 : name="mos_occ_im")
485 66 : alpha = 3.0_dp - REAL(SIZE(mos_old)/2, dp)
486 66 : CALL cp_fm_to_fm(mos_old(re), mos_occ)
487 396 : CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
488 66 : CALL cp_fm_to_fm(mos_old(im), mos_occ_im)
489 396 : CALL cp_fm_column_scale(mos_occ_im, mos(ispin)%occupation_numbers/alpha)
490 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(im)%matrix, &
491 : matrix_v=mos_occ_im, &
492 : matrix_g=mos_occ, &
493 : ncol=ncol, &
494 : alpha=2.0_dp*alpha, &
495 66 : symmetry_mode=-1, keep_sparsity=.FALSE.)
496 :
497 66 : CALL dbcsr_filter(rho_old(im)%matrix, rtp%filter_eps)
498 66 : CALL dbcsr_copy(rho_new(im)%matrix, rho_old(im)%matrix)
499 66 : CALL cp_fm_release(mos_occ_im)
500 172 : CALL cp_fm_release(mos_occ)
501 : END DO
502 : ! Release the mos used to apply the delta kick, no longer required
503 40 : CALL rt_prop_release_mos(rtp)
504 : ELSE
505 50 : DO ispin = 1, SIZE(mos)
506 26 : re = 2*ispin - 1
507 26 : CALL dbcsr_set(rho_old(re)%matrix, 0.0_dp)
508 26 : CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=ncol)
509 :
510 : CALL cp_fm_create(mos_occ, &
511 : matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
512 26 : name="mos_occ")
513 26 : CALL cp_fm_to_fm(mos(ispin)%mo_coeff, mos_occ)
514 26 : IF (mos(ispin)%uniform_occupation) THEN
515 22 : alpha = 3.0_dp - REAL(SIZE(mos), dp)
516 110 : CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
517 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
518 : matrix_v=mos_occ, &
519 : ncol=ncol, &
520 22 : alpha=alpha, keep_sparsity=.FALSE.)
521 : ELSE
522 4 : alpha = 1.0_dp
523 88 : CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
524 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
525 : matrix_v=mos(ispin)%mo_coeff, &
526 : matrix_g=mos_occ, &
527 : ncol=ncol, &
528 4 : alpha=alpha, keep_sparsity=.FALSE.)
529 : END IF
530 26 : CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
531 26 : CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
532 76 : CALL cp_fm_release(mos_occ)
533 : END DO
534 : END IF
535 :
536 64 : END SUBROUTINE rt_initialize_rho_from_mos
537 :
538 : END MODULE rt_propagator_init
|
