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 propagating the orbitals
10 : !> \author Florian Schiffmann (02.09)
11 : ! **************************************************************************************************
12 : MODULE rt_propagation_methods
13 : USE bibliography, ONLY: Kolafa2004,&
14 : Kuhne2007,&
15 : cite_reference
16 : USE cell_types, ONLY: cell_type
17 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_triangular_multiply
18 : USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose
19 : USE cp_cfm_types, ONLY: cp_cfm_create,&
20 : cp_cfm_release,&
21 : cp_cfm_type
22 : USE cp_control_types, ONLY: dft_control_type,&
23 : rtp_control_type
24 : USE cp_dbcsr_api, ONLY: &
25 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
26 : dbcsr_filter, dbcsr_get_block_p, dbcsr_init_p, dbcsr_iterator_blocks_left, &
27 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
28 : dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_transposed, &
29 : dbcsr_type, dbcsr_type_antisymmetric
30 : USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
31 : cp_dbcsr_cholesky_invert
32 : USE cp_dbcsr_contrib, ONLY: dbcsr_frobenius_norm
33 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
34 : dbcsr_allocate_matrix_set,&
35 : dbcsr_deallocate_matrix_set
36 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
37 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
38 : cp_fm_struct_double,&
39 : cp_fm_struct_release,&
40 : cp_fm_struct_type
41 : USE cp_fm_types, ONLY: cp_fm_create,&
42 : cp_fm_get_info,&
43 : cp_fm_release,&
44 : cp_fm_to_fm,&
45 : cp_fm_type
46 : USE cp_log_handling, ONLY: cp_get_default_logger,&
47 : cp_logger_get_default_unit_nr,&
48 : cp_logger_type,&
49 : cp_to_string
50 : USE efield_utils, ONLY: efield_potential_lengh_gauge
51 : USE input_constants, ONLY: do_arnoldi,&
52 : do_bch,&
53 : do_em,&
54 : do_pade,&
55 : do_taylor
56 : USE iterate_matrix, ONLY: matrix_sqrt_Newton_Schulz
57 : USE kinds, ONLY: dp
58 : USE ls_matrix_exp, ONLY: cp_complex_dbcsr_gemm_3
59 : USE mathlib, ONLY: binomial
60 : USE parallel_gemm_api, ONLY: parallel_gemm
61 : USE qs_energy_init, ONLY: qs_energies_init
62 : USE qs_energy_types, ONLY: qs_energy_type
63 : USE qs_environment_types, ONLY: get_qs_env,&
64 : qs_environment_type
65 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
66 : USE qs_ks_types, ONLY: set_ks_env
67 : USE rt_make_propagators, ONLY: propagate_arnoldi,&
68 : propagate_bch,&
69 : propagate_exp,&
70 : propagate_exp_density
71 : USE rt_propagation_output, ONLY: report_density_occupation,&
72 : rt_convergence,&
73 : rt_convergence_density
74 : USE rt_propagation_types, ONLY: get_rtp,&
75 : rt_prop_type
76 : USE rt_propagation_utils, ONLY: calc_S_derivs,&
77 : calc_update_rho,&
78 : calc_update_rho_sparse
79 : USE rt_propagation_velocity_gauge, ONLY: update_vector_potential,&
80 : velocity_gauge_ks_matrix
81 : #include "../base/base_uses.f90"
82 :
83 : IMPLICIT NONE
84 :
85 : PRIVATE
86 :
87 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_methods'
88 :
89 : PUBLIC :: propagation_step, &
90 : s_matrices_create, &
91 : calc_sinvH, &
92 : put_data_to_history
93 :
94 : CONTAINS
95 :
96 : ! **************************************************************************************************
97 : !> \brief performs a single propagation step a(t+Dt)=U(t+Dt,t)*a(0)
98 : !> and calculates the new exponential
99 : !> \param qs_env ...
100 : !> \param rtp ...
101 : !> \param rtp_control ...
102 : !> \author Florian Schiffmann (02.09)
103 : ! **************************************************************************************************
104 :
105 2284 : SUBROUTINE propagation_step(qs_env, rtp, rtp_control)
106 :
107 : TYPE(qs_environment_type), POINTER :: qs_env
108 : TYPE(rt_prop_type), POINTER :: rtp
109 : TYPE(rtp_control_type), POINTER :: rtp_control
110 :
111 : CHARACTER(len=*), PARAMETER :: routineN = 'propagation_step'
112 :
113 : INTEGER :: aspc_order, handle, i, im, re, unit_nr
114 : TYPE(cell_type), POINTER :: cell
115 2284 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: delta_mos, mos_new
116 : TYPE(cp_logger_type), POINTER :: logger
117 2284 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: delta_P, H_last_iter, ks_mix, ks_mix_im, &
118 2284 : matrix_ks, matrix_ks_im, matrix_s, &
119 2284 : rho_new
120 : TYPE(dft_control_type), POINTER :: dft_control
121 :
122 2284 : CALL timeset(routineN, handle)
123 :
124 2284 : logger => cp_get_default_logger()
125 2284 : IF (logger%para_env%is_source()) THEN
126 1142 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
127 : ELSE
128 : unit_nr = -1
129 : END IF
130 :
131 2284 : NULLIFY (cell, delta_P, rho_new, delta_mos, mos_new)
132 2284 : NULLIFY (ks_mix, ks_mix_im)
133 : ! get everything needed and set some values
134 2284 : CALL get_qs_env(qs_env, cell=cell, matrix_s=matrix_s, dft_control=dft_control)
135 :
136 2284 : IF (rtp%iter == 1) THEN
137 614 : CALL qs_energies_init(qs_env, .FALSE.)
138 : !the above recalculates matrix_s, but matrix not changed if ions are fixed
139 614 : IF (rtp_control%fixed_ions) CALL set_ks_env(qs_env%ks_env, s_mstruct_changed=.FALSE.)
140 :
141 : ! add additional terms to matrix_h and matrix_h_im in the case of applied electric field,
142 : ! either in the lengh or velocity gauge.
143 : ! should be called after qs_energies_init and before qs_ks_update_qs_env
144 614 : IF (dft_control%apply_efield_field) THEN
145 216 : IF (ANY(cell%perd(1:3) /= 0)) &
146 0 : CPABORT("Length gauge (efield) and periodicity are not compatible")
147 54 : CALL efield_potential_lengh_gauge(qs_env)
148 560 : ELSE IF (rtp_control%velocity_gauge) THEN
149 32 : IF (dft_control%apply_vector_potential) &
150 32 : CALL update_vector_potential(qs_env, dft_control)
151 32 : CALL velocity_gauge_ks_matrix(qs_env, subtract_nl_term=.FALSE.)
152 : END IF
153 :
154 614 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
155 614 : IF (.NOT. rtp_control%fixed_ions) THEN
156 270 : CALL s_matrices_create(matrix_s, rtp)
157 : END IF
158 614 : rtp%delta_iter = 100.0_dp
159 614 : rtp%mixing_factor = 1.0_dp
160 614 : rtp%mixing = .FALSE.
161 614 : aspc_order = rtp_control%aspc_order
162 614 : CALL aspc_extrapolate(rtp, matrix_s, aspc_order)
163 614 : IF (rtp%linear_scaling) THEN
164 182 : CALL calc_update_rho_sparse(qs_env)
165 : ELSE
166 432 : CALL calc_update_rho(qs_env)
167 : END IF
168 614 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE.)
169 : END IF
170 2284 : IF (.NOT. rtp_control%fixed_ions) THEN
171 1144 : CALL calc_S_derivs(qs_env)
172 : END IF
173 2284 : rtp%converged = .FALSE.
174 :
175 2284 : IF (rtp%linear_scaling) THEN
176 : ! keep temporary copy of the starting density matrix to check for convergence
177 806 : CALL get_rtp(rtp=rtp, rho_new=rho_new)
178 806 : NULLIFY (delta_P)
179 806 : CALL dbcsr_allocate_matrix_set(delta_P, SIZE(rho_new))
180 2910 : DO i = 1, SIZE(rho_new)
181 2104 : CALL dbcsr_init_p(delta_P(i)%matrix)
182 2104 : CALL dbcsr_create(delta_P(i)%matrix, template=rho_new(i)%matrix)
183 2910 : CALL dbcsr_copy(delta_P(i)%matrix, rho_new(i)%matrix)
184 : END DO
185 : ELSE
186 : ! keep temporary copy of the starting mos to check for convergence
187 1478 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
188 8214 : ALLOCATE (delta_mos(SIZE(mos_new)))
189 5258 : DO i = 1, SIZE(mos_new)
190 : CALL cp_fm_create(delta_mos(i), &
191 : matrix_struct=mos_new(i)%matrix_struct, &
192 3780 : name="delta_mos"//TRIM(ADJUSTL(cp_to_string(i))))
193 5258 : CALL cp_fm_to_fm(mos_new(i), delta_mos(i))
194 : END DO
195 : END IF
196 :
197 : CALL get_qs_env(qs_env, &
198 : matrix_ks=matrix_ks, &
199 2284 : matrix_ks_im=matrix_ks_im)
200 :
201 2284 : CALL get_rtp(rtp=rtp, H_last_iter=H_last_iter)
202 2284 : IF (rtp%mixing) THEN
203 96 : IF (unit_nr > 0) THEN
204 48 : WRITE (unit_nr, '(t3,a,2f16.8)') "Mixing the Hamiltonians to improve robustness, mixing factor: ", rtp%mixing_factor
205 : END IF
206 96 : CALL dbcsr_allocate_matrix_set(ks_mix, SIZE(matrix_ks))
207 96 : CALL dbcsr_allocate_matrix_set(ks_mix_im, SIZE(matrix_ks))
208 192 : DO i = 1, SIZE(matrix_ks)
209 96 : CALL dbcsr_init_p(ks_mix(i)%matrix)
210 96 : CALL dbcsr_create(ks_mix(i)%matrix, template=matrix_ks(1)%matrix)
211 96 : CALL dbcsr_init_p(ks_mix_im(i)%matrix)
212 192 : CALL dbcsr_create(ks_mix_im(i)%matrix, template=matrix_ks(1)%matrix, matrix_type=dbcsr_type_antisymmetric)
213 : END DO
214 192 : DO i = 1, SIZE(matrix_ks)
215 96 : re = 2*i - 1
216 96 : im = 2*i
217 96 : CALL dbcsr_add(ks_mix(i)%matrix, matrix_ks(i)%matrix, 0.0_dp, rtp%mixing_factor)
218 96 : CALL dbcsr_add(ks_mix(i)%matrix, H_last_iter(re)%matrix, 1.0_dp, 1.0_dp - rtp%mixing_factor)
219 192 : IF (rtp%propagate_complex_ks) THEN
220 0 : CALL dbcsr_add(ks_mix_im(i)%matrix, matrix_ks_im(i)%matrix, 0.0_dp, rtp%mixing_factor)
221 0 : CALL dbcsr_add(ks_mix_im(i)%matrix, H_last_iter(im)%matrix, 1.0_dp, 1.0_dp - rtp%mixing_factor)
222 : END IF
223 : END DO
224 96 : CALL calc_SinvH(rtp, ks_mix, ks_mix_im, rtp_control)
225 192 : DO i = 1, SIZE(matrix_ks)
226 96 : re = 2*i - 1
227 96 : im = 2*i
228 96 : CALL dbcsr_copy(H_last_iter(re)%matrix, ks_mix(i)%matrix)
229 192 : IF (rtp%propagate_complex_ks) THEN
230 0 : CALL dbcsr_copy(H_last_iter(im)%matrix, ks_mix_im(i)%matrix)
231 : END IF
232 : END DO
233 96 : CALL dbcsr_deallocate_matrix_set(ks_mix)
234 96 : CALL dbcsr_deallocate_matrix_set(ks_mix_im)
235 : ELSE
236 2188 : CALL calc_SinvH(rtp, matrix_ks, matrix_ks_im, rtp_control)
237 5034 : DO i = 1, SIZE(matrix_ks)
238 2846 : re = 2*i - 1
239 2846 : im = 2*i
240 2846 : CALL dbcsr_copy(H_last_iter(re)%matrix, matrix_ks(i)%matrix)
241 5034 : IF (rtp%propagate_complex_ks) THEN
242 422 : CALL dbcsr_copy(H_last_iter(im)%matrix, matrix_ks_im(i)%matrix)
243 : END IF
244 : END DO
245 : END IF
246 :
247 2284 : CALL compute_propagator_matrix(rtp, rtp_control%propagator)
248 :
249 3960 : SELECT CASE (rtp_control%mat_exp)
250 : CASE (do_pade, do_taylor)
251 1676 : IF (rtp%linear_scaling) THEN
252 622 : CALL propagate_exp_density(rtp, rtp_control)
253 622 : CALL calc_update_rho_sparse(qs_env)
254 : ELSE
255 1054 : CALL propagate_exp(rtp, rtp_control)
256 1054 : CALL calc_update_rho(qs_env)
257 : END IF
258 : CASE (do_arnoldi)
259 424 : CALL propagate_arnoldi(rtp, rtp_control)
260 424 : CALL calc_update_rho(qs_env)
261 : CASE (do_bch)
262 184 : CALL propagate_bch(rtp, rtp_control)
263 2468 : CALL calc_update_rho_sparse(qs_env)
264 : END SELECT
265 2284 : CALL step_finalize(qs_env, rtp_control, delta_mos, delta_P)
266 2284 : IF (rtp%linear_scaling) THEN
267 806 : CALL dbcsr_deallocate_matrix_set(delta_P)
268 : ELSE
269 1478 : CALL cp_fm_release(delta_mos)
270 : END IF
271 :
272 2284 : CALL timestop(handle)
273 :
274 2284 : END SUBROUTINE propagation_step
275 :
276 : ! **************************************************************************************************
277 : !> \brief Performs all the stuff to finish the step:
278 : !> convergence checks
279 : !> copying stuff into right place for the next step
280 : !> updating the history for extrapolation
281 : !> \param qs_env ...
282 : !> \param rtp_control ...
283 : !> \param delta_mos ...
284 : !> \param delta_P ...
285 : !> \author Florian Schiffmann (02.09)
286 : ! **************************************************************************************************
287 :
288 2284 : SUBROUTINE step_finalize(qs_env, rtp_control, delta_mos, delta_P)
289 : TYPE(qs_environment_type), POINTER :: qs_env
290 : TYPE(rtp_control_type), POINTER :: rtp_control
291 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: delta_mos
292 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: delta_P
293 :
294 : CHARACTER(len=*), PARAMETER :: routineN = 'step_finalize'
295 :
296 : INTEGER :: handle, i, ihist
297 2284 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new, mos_old
298 2284 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: exp_H_new, exp_H_old, matrix_ks, &
299 2284 : matrix_ks_im, rho_new, rho_old, s_mat
300 : TYPE(qs_energy_type), POINTER :: energy
301 : TYPE(rt_prop_type), POINTER :: rtp
302 :
303 2284 : CALL timeset(routineN, handle)
304 :
305 : CALL get_qs_env(qs_env=qs_env, rtp=rtp, matrix_s=s_mat, &
306 2284 : matrix_ks=matrix_ks, matrix_ks_im=matrix_ks_im, energy=energy)
307 2284 : CALL get_rtp(rtp=rtp, exp_H_old=exp_H_old, exp_H_new=exp_H_new)
308 :
309 2284 : IF (rtp_control%sc_check_start .LT. rtp%iter) THEN
310 2284 : rtp%delta_iter_old = rtp%delta_iter
311 2284 : IF (rtp%linear_scaling) THEN
312 806 : CALL rt_convergence_density(rtp, delta_P, rtp%delta_iter)
313 : ELSE
314 1478 : CALL rt_convergence(rtp, s_mat(1)%matrix, delta_mos, rtp%delta_iter)
315 : END IF
316 2284 : rtp%converged = (rtp%delta_iter .LT. rtp_control%eps_ener)
317 : !Apply mixing if scf loop is not converging
318 :
319 : !It would be better to redo the current step with mixixng,
320 : !but currently the decision is made to use mixing from the next step on
321 2284 : IF (rtp_control%sc_check_start .LT. rtp%iter + 1) THEN
322 2284 : IF (rtp%delta_iter/rtp%delta_iter_old > 0.9) THEN
323 6 : rtp%mixing_factor = MAX(rtp%mixing_factor/2.0_dp, 0.125_dp)
324 6 : rtp%mixing = .TRUE.
325 : END IF
326 : END IF
327 : END IF
328 :
329 2284 : IF (rtp%converged) THEN
330 614 : IF (rtp%linear_scaling) THEN
331 182 : CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
332 : CALL purify_mcweeny_complex_nonorth(rho_new, s_mat, rtp%filter_eps, rtp%filter_eps_small, &
333 182 : rtp_control%mcweeny_max_iter, rtp_control%mcweeny_eps)
334 182 : IF (rtp_control%mcweeny_max_iter > 0) CALL calc_update_rho_sparse(qs_env)
335 182 : CALL report_density_occupation(rtp%filter_eps, rho_new)
336 686 : DO i = 1, SIZE(rho_new)
337 686 : CALL dbcsr_copy(rho_old(i)%matrix, rho_new(i)%matrix)
338 : END DO
339 : ELSE
340 432 : CALL get_rtp(rtp=rtp, mos_old=mos_old, mos_new=mos_new)
341 1500 : DO i = 1, SIZE(mos_new)
342 1500 : CALL cp_fm_to_fm(mos_new(i), mos_old(i))
343 : END DO
344 : END IF
345 614 : IF (rtp_control%propagator == do_em) CALL calc_SinvH(rtp, matrix_ks, matrix_ks_im, rtp_control)
346 2186 : DO i = 1, SIZE(exp_H_new)
347 2186 : CALL dbcsr_copy(exp_H_old(i)%matrix, exp_H_new(i)%matrix)
348 : END DO
349 614 : ihist = MOD(rtp%istep, rtp_control%aspc_order) + 1
350 614 : IF (rtp_control%fixed_ions) THEN
351 344 : CALL put_data_to_history(rtp, rho=rho_new, mos=mos_new, ihist=ihist)
352 : ELSE
353 270 : CALL put_data_to_history(rtp, rho=rho_new, mos=mos_new, s_mat=s_mat, ihist=ihist)
354 : END IF
355 : END IF
356 :
357 2284 : rtp%energy_new = energy%total
358 :
359 2284 : CALL timestop(handle)
360 :
361 2284 : END SUBROUTINE step_finalize
362 :
363 : ! **************************************************************************************************
364 : !> \brief computes the propagator matrix for EM/ETRS, RTP/EMD
365 : !> \param rtp ...
366 : !> \param propagator ...
367 : !> \author Florian Schiffmann (02.09)
368 : ! **************************************************************************************************
369 :
370 4568 : SUBROUTINE compute_propagator_matrix(rtp, propagator)
371 : TYPE(rt_prop_type), POINTER :: rtp
372 : INTEGER :: propagator
373 :
374 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_propagator_matrix'
375 :
376 : INTEGER :: handle, i
377 : REAL(Kind=dp) :: dt, prefac
378 2284 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: exp_H_new, exp_H_old, propagator_matrix
379 :
380 2284 : CALL timeset(routineN, handle)
381 : CALL get_rtp(rtp=rtp, exp_H_new=exp_H_new, exp_H_old=exp_H_old, &
382 2284 : propagator_matrix=propagator_matrix, dt=dt)
383 :
384 2284 : prefac = -0.5_dp*dt
385 :
386 8168 : DO i = 1, SIZE(exp_H_new)
387 5884 : CALL dbcsr_add(propagator_matrix(i)%matrix, exp_H_new(i)%matrix, 0.0_dp, prefac)
388 5884 : IF (propagator == do_em) &
389 2348 : CALL dbcsr_add(propagator_matrix(i)%matrix, exp_H_old(i)%matrix, 1.0_dp, prefac)
390 : END DO
391 :
392 2284 : CALL timestop(handle)
393 :
394 2284 : END SUBROUTINE compute_propagator_matrix
395 :
396 : ! **************************************************************************************************
397 : !> \brief computes S_inv*H, if needed Sinv*B and S_inv*H_imag and store these quantities to the
398 : !> \brief exp_H for the real and imag part (for RTP and EMD)
399 : !> \param rtp ...
400 : !> \param matrix_ks ...
401 : !> \param matrix_ks_im ...
402 : !> \param rtp_control ...
403 : !> \author Florian Schiffmann (02.09)
404 : ! **************************************************************************************************
405 :
406 2494 : SUBROUTINE calc_SinvH(rtp, matrix_ks, matrix_ks_im, rtp_control)
407 : TYPE(rt_prop_type), POINTER :: rtp
408 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_im
409 : TYPE(rtp_control_type), POINTER :: rtp_control
410 :
411 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_SinvH'
412 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
413 :
414 : INTEGER :: handle, im, ispin, re
415 2494 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: exp_H, SinvB, SinvH, SinvH_imag
416 : TYPE(dbcsr_type) :: matrix_ks_nosym
417 : TYPE(dbcsr_type), POINTER :: B_mat, S_inv
418 :
419 2494 : CALL timeset(routineN, handle)
420 2494 : CALL get_rtp(rtp=rtp, S_inv=S_inv, exp_H_new=exp_H)
421 5714 : DO ispin = 1, SIZE(matrix_ks)
422 3220 : re = ispin*2 - 1
423 3220 : im = ispin*2
424 3220 : CALL dbcsr_set(exp_H(re)%matrix, zero)
425 5714 : CALL dbcsr_set(exp_H(im)%matrix, zero)
426 : END DO
427 2494 : CALL dbcsr_create(matrix_ks_nosym, template=matrix_ks(1)%matrix, matrix_type="N")
428 :
429 : ! Real part of S_inv x H -> imag part of exp_H
430 5714 : DO ispin = 1, SIZE(matrix_ks)
431 3220 : re = ispin*2 - 1
432 3220 : im = ispin*2
433 3220 : CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, matrix_ks_nosym)
434 : CALL dbcsr_multiply("N", "N", one, S_inv, matrix_ks_nosym, zero, exp_H(im)%matrix, &
435 3220 : filter_eps=rtp%filter_eps)
436 5714 : IF (.NOT. rtp_control%fixed_ions) THEN
437 1466 : CALL get_rtp(rtp=rtp, SinvH=SinvH)
438 1466 : CALL dbcsr_copy(SinvH(ispin)%matrix, exp_H(im)%matrix)
439 : END IF
440 : END DO
441 :
442 : ! Imag part of S_inv x H -> real part of exp_H
443 2494 : IF (rtp%propagate_complex_ks) THEN
444 910 : DO ispin = 1, SIZE(matrix_ks)
445 466 : re = ispin*2 - 1
446 466 : im = ispin*2
447 466 : CALL dbcsr_set(matrix_ks_nosym, 0.0_dp)
448 466 : CALL dbcsr_desymmetrize(matrix_ks_im(ispin)%matrix, matrix_ks_nosym)
449 : ! - SinvH_imag is added to exp_H(re)%matrix
450 : CALL dbcsr_multiply("N", "N", -one, S_inv, matrix_ks_nosym, zero, exp_H(re)%matrix, &
451 466 : filter_eps=rtp%filter_eps)
452 910 : IF (.NOT. rtp_control%fixed_ions) THEN
453 266 : CALL get_rtp(rtp=rtp, SinvH_imag=SinvH_imag)
454 : ! -SinvH_imag is saved
455 266 : CALL dbcsr_copy(SinvH_imag(ispin)%matrix, exp_H(re)%matrix)
456 : END IF
457 : END DO
458 : END IF
459 : ! EMD case: the real part of exp_H should be updated with B
460 2494 : IF (.NOT. rtp_control%fixed_ions) THEN
461 1224 : CALL get_rtp(rtp=rtp, B_mat=B_mat, SinvB=SinvB)
462 1224 : CALL dbcsr_set(matrix_ks_nosym, 0.0_dp)
463 1224 : CALL dbcsr_multiply("N", "N", one, S_inv, B_mat, zero, matrix_ks_nosym, filter_eps=rtp%filter_eps)
464 2690 : DO ispin = 1, SIZE(matrix_ks)
465 1466 : re = ispin*2 - 1
466 1466 : im = ispin*2
467 : ! + SinvB is added to exp_H(re)%matrix
468 1466 : CALL dbcsr_add(exp_H(re)%matrix, matrix_ks_nosym, 1.0_dp, 1.0_dp)
469 : ! + SinvB is saved
470 2690 : CALL dbcsr_copy(SinvB(ispin)%matrix, matrix_ks_nosym)
471 : END DO
472 : END IF
473 : ! Otherwise no real part for exp_H
474 :
475 2494 : CALL dbcsr_release(matrix_ks_nosym)
476 2494 : CALL timestop(handle)
477 :
478 2494 : END SUBROUTINE calc_SinvH
479 :
480 : ! **************************************************************************************************
481 : !> \brief calculates the needed overlap-like matrices
482 : !> depending on the way the exponential is calculated, only S^-1 is needed
483 : !> \param s_mat ...
484 : !> \param rtp ...
485 : !> \author Florian Schiffmann (02.09)
486 : ! **************************************************************************************************
487 :
488 468 : SUBROUTINE s_matrices_create(s_mat, rtp)
489 :
490 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: s_mat
491 : TYPE(rt_prop_type), POINTER :: rtp
492 :
493 : CHARACTER(len=*), PARAMETER :: routineN = 's_matrices_create'
494 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
495 :
496 : INTEGER :: handle
497 : TYPE(dbcsr_type), POINTER :: S_half, S_inv, S_minus_half
498 :
499 468 : CALL timeset(routineN, handle)
500 :
501 468 : CALL get_rtp(rtp=rtp, S_inv=S_inv)
502 :
503 468 : IF (rtp%linear_scaling) THEN
504 136 : CALL get_rtp(rtp=rtp, S_half=S_half, S_minus_half=S_minus_half)
505 : CALL matrix_sqrt_Newton_Schulz(S_half, S_minus_half, s_mat(1)%matrix, rtp%filter_eps, &
506 136 : rtp%newton_schulz_order, rtp%lanzcos_threshold, rtp%lanzcos_max_iter)
507 : CALL dbcsr_multiply("N", "N", one, S_minus_half, S_minus_half, zero, S_inv, &
508 136 : filter_eps=rtp%filter_eps)
509 : ELSE
510 332 : CALL dbcsr_copy(S_inv, s_mat(1)%matrix)
511 : CALL cp_dbcsr_cholesky_decompose(S_inv, para_env=rtp%ao_ao_fmstruct%para_env, &
512 332 : blacs_env=rtp%ao_ao_fmstruct%context)
513 : CALL cp_dbcsr_cholesky_invert(S_inv, para_env=rtp%ao_ao_fmstruct%para_env, &
514 332 : blacs_env=rtp%ao_ao_fmstruct%context, uplo_to_full=.TRUE.)
515 : END IF
516 :
517 468 : CALL timestop(handle)
518 468 : END SUBROUTINE s_matrices_create
519 :
520 : ! **************************************************************************************************
521 : !> \brief Calculates the frobenius norm of a complex matrix represented by two real matrices
522 : !> \param frob_norm ...
523 : !> \param mat_re ...
524 : !> \param mat_im ...
525 : !> \author Samuel Andermatt (04.14)
526 : ! **************************************************************************************************
527 :
528 568 : SUBROUTINE complex_frobenius_norm(frob_norm, mat_re, mat_im)
529 :
530 : REAL(KIND=dp), INTENT(out) :: frob_norm
531 : TYPE(dbcsr_type), POINTER :: mat_re, mat_im
532 :
533 : CHARACTER(len=*), PARAMETER :: routineN = 'complex_frobenius_norm'
534 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
535 :
536 : INTEGER :: col_atom, handle, row_atom
537 : LOGICAL :: found
538 284 : REAL(dp), DIMENSION(:, :), POINTER :: block_values, block_values2
539 : TYPE(dbcsr_iterator_type) :: iter
540 : TYPE(dbcsr_type), POINTER :: tmp
541 :
542 284 : CALL timeset(routineN, handle)
543 :
544 : NULLIFY (tmp)
545 284 : ALLOCATE (tmp)
546 284 : CALL dbcsr_create(tmp, template=mat_re)
547 : !make sure the tmp has the same sparsity pattern as the real and the complex part combined
548 284 : CALL dbcsr_add(tmp, mat_re, zero, one)
549 284 : CALL dbcsr_add(tmp, mat_im, zero, one)
550 284 : CALL dbcsr_set(tmp, zero)
551 : !calculate the hadamard product
552 284 : CALL dbcsr_iterator_start(iter, tmp)
553 1804 : DO WHILE (dbcsr_iterator_blocks_left(iter))
554 1520 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
555 1520 : CALL dbcsr_get_block_p(mat_re, row_atom, col_atom, block_values2, found=found)
556 1520 : IF (found) THEN
557 264788 : block_values = block_values2*block_values2
558 : END IF
559 1520 : CALL dbcsr_get_block_p(mat_im, row_atom, col_atom, block_values2, found=found)
560 1520 : IF (found) THEN
561 250308 : block_values = block_values + block_values2*block_values2
562 : END IF
563 133438 : block_values = SQRT(block_values)
564 : END DO
565 284 : CALL dbcsr_iterator_stop(iter)
566 284 : frob_norm = dbcsr_frobenius_norm(tmp)
567 :
568 284 : CALL dbcsr_deallocate_matrix(tmp)
569 :
570 284 : CALL timestop(handle)
571 :
572 284 : END SUBROUTINE complex_frobenius_norm
573 :
574 : ! **************************************************************************************************
575 : !> \brief Does McWeeny for complex matrices in the non-orthogonal basis
576 : !> \param P ...
577 : !> \param s_mat ...
578 : !> \param eps ...
579 : !> \param eps_small ...
580 : !> \param max_iter ...
581 : !> \param threshold ...
582 : !> \author Samuel Andermatt (04.14)
583 : ! **************************************************************************************************
584 :
585 182 : SUBROUTINE purify_mcweeny_complex_nonorth(P, s_mat, eps, eps_small, max_iter, threshold)
586 :
587 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: P, s_mat
588 : REAL(KIND=dp), INTENT(in) :: eps, eps_small
589 : INTEGER, INTENT(in) :: max_iter
590 : REAL(KIND=dp), INTENT(in) :: threshold
591 :
592 : CHARACTER(len=*), PARAMETER :: routineN = 'purify_mcweeny_complex_nonorth'
593 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
594 :
595 : INTEGER :: handle, i, im, imax, ispin, re, unit_nr
596 : REAL(KIND=dp) :: frob_norm
597 : TYPE(cp_logger_type), POINTER :: logger
598 182 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: PS, PSP, tmp
599 :
600 182 : CALL timeset(routineN, handle)
601 :
602 182 : logger => cp_get_default_logger()
603 182 : IF (logger%para_env%is_source()) THEN
604 91 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
605 : ELSE
606 : unit_nr = -1
607 : END IF
608 :
609 182 : NULLIFY (tmp, PS, PSP)
610 182 : CALL dbcsr_allocate_matrix_set(tmp, SIZE(P))
611 182 : CALL dbcsr_allocate_matrix_set(PSP, SIZE(P))
612 182 : CALL dbcsr_allocate_matrix_set(PS, SIZE(P))
613 686 : DO i = 1, SIZE(P)
614 504 : CALL dbcsr_init_p(PS(i)%matrix)
615 504 : CALL dbcsr_create(PS(i)%matrix, template=P(1)%matrix)
616 504 : CALL dbcsr_init_p(PSP(i)%matrix)
617 504 : CALL dbcsr_create(PSP(i)%matrix, template=P(1)%matrix)
618 504 : CALL dbcsr_init_p(tmp(i)%matrix)
619 686 : CALL dbcsr_create(tmp(i)%matrix, template=P(1)%matrix)
620 : END DO
621 182 : IF (SIZE(P) == 2) THEN
622 112 : CALL dbcsr_scale(P(1)%matrix, one/2)
623 112 : CALL dbcsr_scale(P(2)%matrix, one/2)
624 : END IF
625 434 : DO ispin = 1, SIZE(P)/2
626 252 : re = 2*ispin - 1
627 252 : im = 2*ispin
628 252 : imax = MAX(max_iter, 1) !if max_iter is 0 then only the deviation from idempotency needs to be calculated
629 516 : DO i = 1, imax
630 : CALL dbcsr_multiply("N", "N", one, P(re)%matrix, s_mat(1)%matrix, &
631 284 : zero, PS(re)%matrix, filter_eps=eps_small)
632 : CALL dbcsr_multiply("N", "N", one, P(im)%matrix, s_mat(1)%matrix, &
633 284 : zero, PS(im)%matrix, filter_eps=eps_small)
634 : CALL cp_complex_dbcsr_gemm_3("N", "N", one, PS(re)%matrix, PS(im)%matrix, &
635 : P(re)%matrix, P(im)%matrix, zero, PSP(re)%matrix, PSP(im)%matrix, &
636 284 : filter_eps=eps_small)
637 284 : CALL dbcsr_copy(tmp(re)%matrix, PSP(re)%matrix)
638 284 : CALL dbcsr_copy(tmp(im)%matrix, PSP(im)%matrix)
639 284 : CALL dbcsr_add(tmp(re)%matrix, P(re)%matrix, 1.0_dp, -1.0_dp)
640 284 : CALL dbcsr_add(tmp(im)%matrix, P(im)%matrix, 1.0_dp, -1.0_dp)
641 284 : CALL complex_frobenius_norm(frob_norm, tmp(re)%matrix, tmp(im)%matrix)
642 284 : IF (unit_nr .GT. 0) WRITE (unit_nr, '(t3,a,2f16.8)') "Deviation from idempotency: ", frob_norm
643 800 : IF (frob_norm .GT. threshold .AND. max_iter > 0) THEN
644 264 : CALL dbcsr_copy(P(re)%matrix, PSP(re)%matrix)
645 264 : CALL dbcsr_copy(P(im)%matrix, PSP(im)%matrix)
646 : CALL cp_complex_dbcsr_gemm_3("N", "N", -2.0_dp, PS(re)%matrix, PS(im)%matrix, &
647 : PSP(re)%matrix, PSP(im)%matrix, 3.0_dp, P(re)%matrix, P(im)%matrix, &
648 264 : filter_eps=eps_small)
649 264 : CALL dbcsr_filter(P(re)%matrix, eps)
650 264 : CALL dbcsr_filter(P(im)%matrix, eps)
651 : !make sure P is exactly hermitian
652 264 : CALL dbcsr_transposed(tmp(re)%matrix, P(re)%matrix)
653 264 : CALL dbcsr_add(P(re)%matrix, tmp(re)%matrix, one/2, one/2)
654 264 : CALL dbcsr_transposed(tmp(im)%matrix, P(im)%matrix)
655 264 : CALL dbcsr_add(P(im)%matrix, tmp(im)%matrix, one/2, -one/2)
656 : ELSE
657 : EXIT
658 : END IF
659 : END DO
660 : !make sure P is hermitian
661 252 : CALL dbcsr_transposed(tmp(re)%matrix, P(re)%matrix)
662 252 : CALL dbcsr_add(P(re)%matrix, tmp(re)%matrix, one/2, one/2)
663 252 : CALL dbcsr_transposed(tmp(im)%matrix, P(im)%matrix)
664 434 : CALL dbcsr_add(P(im)%matrix, tmp(im)%matrix, one/2, -one/2)
665 : END DO
666 182 : IF (SIZE(P) == 2) THEN
667 112 : CALL dbcsr_scale(P(1)%matrix, one*2)
668 112 : CALL dbcsr_scale(P(2)%matrix, one*2)
669 : END IF
670 182 : CALL dbcsr_deallocate_matrix_set(tmp)
671 182 : CALL dbcsr_deallocate_matrix_set(PS)
672 182 : CALL dbcsr_deallocate_matrix_set(PSP)
673 :
674 182 : CALL timestop(handle)
675 :
676 182 : END SUBROUTINE purify_mcweeny_complex_nonorth
677 :
678 : ! **************************************************************************************************
679 : !> \brief ...
680 : !> \param rtp ...
681 : !> \param matrix_s ...
682 : !> \param aspc_order ...
683 : ! **************************************************************************************************
684 614 : SUBROUTINE aspc_extrapolate(rtp, matrix_s, aspc_order)
685 : TYPE(rt_prop_type), POINTER :: rtp
686 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
687 : INTEGER, INTENT(in) :: aspc_order
688 :
689 : CHARACTER(len=*), PARAMETER :: routineN = 'aspc_extrapolate'
690 : COMPLEX(KIND=dp), PARAMETER :: cone = (1.0_dp, 0.0_dp), &
691 : czero = (0.0_dp, 0.0_dp)
692 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
693 :
694 : INTEGER :: handle, i, iaspc, icol_local, ihist, &
695 : imat, k, kdbl, n, naspc, ncol_local, &
696 : nmat
697 : REAL(KIND=dp) :: alpha
698 : TYPE(cp_cfm_type) :: cfm_tmp, cfm_tmp1, csc
699 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct, matrix_struct_new
700 : TYPE(cp_fm_type) :: fm_tmp, fm_tmp1, fm_tmp2
701 614 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
702 614 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: mo_hist
703 614 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new, s_hist
704 614 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_hist
705 :
706 614 : NULLIFY (rho_hist)
707 614 : CALL timeset(routineN, handle)
708 614 : CALL cite_reference(Kolafa2004)
709 614 : CALL cite_reference(Kuhne2007)
710 :
711 614 : IF (rtp%linear_scaling) THEN
712 182 : CALL get_rtp(rtp=rtp, rho_new=rho_new)
713 : ELSE
714 432 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
715 : END IF
716 :
717 614 : naspc = MIN(rtp%istep, aspc_order)
718 614 : IF (rtp%linear_scaling) THEN
719 182 : nmat = SIZE(rho_new)
720 182 : rho_hist => rtp%history%rho_history
721 686 : DO imat = 1, nmat
722 1454 : DO iaspc = 1, naspc
723 : alpha = (-1.0_dp)**(iaspc + 1)*REAL(iaspc, KIND=dp)* &
724 768 : binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1)
725 768 : ihist = MOD(rtp%istep - iaspc, aspc_order) + 1
726 1272 : IF (iaspc == 1) THEN
727 504 : CALL dbcsr_add(rho_new(imat)%matrix, rho_hist(imat, ihist)%matrix, zero, alpha)
728 : ELSE
729 264 : CALL dbcsr_add(rho_new(imat)%matrix, rho_hist(imat, ihist)%matrix, one, alpha)
730 : END IF
731 : END DO
732 : END DO
733 : ELSE
734 432 : mo_hist => rtp%history%mo_history
735 432 : nmat = SIZE(mos_new)
736 1500 : DO imat = 1, nmat
737 3876 : DO iaspc = 1, naspc
738 : alpha = (-1.0_dp)**(iaspc + 1)*REAL(iaspc, KIND=dp)* &
739 2376 : binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1)
740 2376 : ihist = MOD(rtp%istep - iaspc, aspc_order) + 1
741 3444 : IF (iaspc == 1) THEN
742 1068 : CALL cp_fm_scale_and_add(zero, mos_new(imat), alpha, mo_hist(imat, ihist))
743 : ELSE
744 1308 : CALL cp_fm_scale_and_add(one, mos_new(imat), alpha, mo_hist(imat, ihist))
745 : END IF
746 : END DO
747 : END DO
748 :
749 432 : mo_hist => rtp%history%mo_history
750 432 : s_hist => rtp%history%s_history
751 966 : DO i = 1, SIZE(mos_new)/2
752 534 : NULLIFY (matrix_struct, matrix_struct_new)
753 :
754 : CALL cp_fm_struct_double(matrix_struct, &
755 : mos_new(2*i)%matrix_struct, &
756 : mos_new(2*i)%matrix_struct%context, &
757 534 : .TRUE., .FALSE.)
758 :
759 534 : CALL cp_fm_create(fm_tmp, matrix_struct)
760 534 : CALL cp_fm_create(fm_tmp1, matrix_struct)
761 534 : CALL cp_fm_create(fm_tmp2, mos_new(2*i)%matrix_struct)
762 534 : CALL cp_cfm_create(cfm_tmp, mos_new(2*i)%matrix_struct)
763 534 : CALL cp_cfm_create(cfm_tmp1, mos_new(2*i)%matrix_struct)
764 :
765 534 : CALL cp_fm_get_info(fm_tmp, ncol_global=kdbl)
766 :
767 : CALL cp_fm_get_info(mos_new(2*i), &
768 : nrow_global=n, &
769 : ncol_global=k, &
770 534 : ncol_local=ncol_local)
771 :
772 : CALL cp_fm_struct_create(matrix_struct_new, &
773 : template_fmstruct=mos_new(2*i)%matrix_struct, &
774 : nrow_global=k, &
775 534 : ncol_global=k)
776 534 : CALL cp_cfm_create(csc, matrix_struct_new)
777 :
778 534 : CALL cp_fm_struct_release(matrix_struct_new)
779 534 : CALL cp_fm_struct_release(matrix_struct)
780 :
781 : ! first the most recent
782 :
783 : ! reorthogonalize vectors
784 2292 : DO icol_local = 1, ncol_local
785 13980 : fm_tmp%local_data(:, icol_local) = mos_new(2*i - 1)%local_data(:, icol_local)
786 14514 : fm_tmp%local_data(:, icol_local + ncol_local) = mos_new(2*i)%local_data(:, icol_local)
787 : END DO
788 :
789 534 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, fm_tmp, fm_tmp1, kdbl)
790 :
791 2292 : DO icol_local = 1, ncol_local
792 : cfm_tmp%local_data(:, icol_local) = CMPLX(fm_tmp1%local_data(:, icol_local), &
793 13980 : fm_tmp1%local_data(:, icol_local + ncol_local), dp)
794 : cfm_tmp1%local_data(:, icol_local) = CMPLX(mos_new(2*i - 1)%local_data(:, icol_local), &
795 14514 : mos_new(2*i)%local_data(:, icol_local), dp)
796 : END DO
797 534 : CALL parallel_gemm('C', 'N', k, k, n, cone, cfm_tmp1, cfm_tmp, czero, csc)
798 534 : CALL cp_cfm_cholesky_decompose(csc)
799 534 : CALL cp_cfm_triangular_multiply(csc, cfm_tmp1, n_cols=k, side='R', invert_tr=.TRUE.)
800 2292 : DO icol_local = 1, ncol_local
801 13980 : mos_new(2*i - 1)%local_data(:, icol_local) = REAL(cfm_tmp1%local_data(:, icol_local), dp)
802 14514 : mos_new(2*i)%local_data(:, icol_local) = AIMAG(cfm_tmp1%local_data(:, icol_local))
803 : END DO
804 :
805 : ! deallocate work matrices
806 534 : CALL cp_cfm_release(csc)
807 534 : CALL cp_fm_release(fm_tmp)
808 534 : CALL cp_fm_release(fm_tmp1)
809 534 : CALL cp_fm_release(fm_tmp2)
810 534 : CALL cp_cfm_release(cfm_tmp)
811 2034 : CALL cp_cfm_release(cfm_tmp1)
812 : END DO
813 :
814 : END IF
815 :
816 614 : CALL timestop(handle)
817 :
818 614 : END SUBROUTINE aspc_extrapolate
819 :
820 : ! **************************************************************************************************
821 : !> \brief ...
822 : !> \param rtp ...
823 : !> \param mos ...
824 : !> \param rho ...
825 : !> \param s_mat ...
826 : !> \param ihist ...
827 : ! **************************************************************************************************
828 812 : SUBROUTINE put_data_to_history(rtp, mos, rho, s_mat, ihist)
829 : TYPE(rt_prop_type), POINTER :: rtp
830 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos
831 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho
832 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
833 : POINTER :: s_mat
834 : INTEGER :: ihist
835 :
836 : INTEGER :: i
837 :
838 812 : IF (rtp%linear_scaling) THEN
839 1032 : DO i = 1, SIZE(rho)
840 1032 : CALL dbcsr_copy(rtp%history%rho_history(i, ihist)%matrix, rho(i)%matrix)
841 : END DO
842 : ELSE
843 1884 : DO i = 1, SIZE(mos)
844 1884 : CALL cp_fm_to_fm(mos(i), rtp%history%mo_history(i, ihist))
845 : END DO
846 540 : IF (PRESENT(s_mat)) THEN
847 332 : IF (ASSOCIATED(rtp%history%s_history(ihist)%matrix)) THEN ! the sparsity might be different
848 : ! (future struct:check)
849 124 : CALL dbcsr_deallocate_matrix(rtp%history%s_history(ihist)%matrix)
850 : END IF
851 332 : ALLOCATE (rtp%history%s_history(ihist)%matrix)
852 332 : CALL dbcsr_copy(rtp%history%s_history(ihist)%matrix, s_mat(1)%matrix)
853 : END IF
854 : END IF
855 :
856 812 : END SUBROUTINE put_data_to_history
857 :
858 : END MODULE rt_propagation_methods
|
