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