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