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