Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 : Schreder2021,&
16 : cite_reference
17 : USE cell_types, ONLY: cell_type
18 : USE cp_array_utils, ONLY: cp_1d_r_p_type
19 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_triangular_multiply
20 : USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose
21 : USE cp_cfm_types, ONLY: cp_cfm_create,&
22 : cp_cfm_release,&
23 : cp_cfm_type
24 : USE cp_control_types, ONLY: dft_control_type,&
25 : rtp_control_type
26 : USE cp_dbcsr_api, ONLY: &
27 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
28 : dbcsr_filter, dbcsr_get_block_p, dbcsr_init_p, dbcsr_iterator_blocks_left, &
29 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
30 : dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_transposed, &
31 : dbcsr_type, dbcsr_type_antisymmetric
32 : USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
33 : cp_dbcsr_cholesky_invert
34 : USE cp_dbcsr_contrib, ONLY: dbcsr_frobenius_norm
35 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
36 : dbcsr_allocate_matrix_set,&
37 : dbcsr_deallocate_matrix_set
38 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
39 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
40 : cp_fm_struct_double,&
41 : cp_fm_struct_release,&
42 : cp_fm_struct_type
43 : USE cp_fm_types, ONLY: cp_fm_create,&
44 : cp_fm_get_info,&
45 : cp_fm_release,&
46 : cp_fm_to_fm,&
47 : cp_fm_type
48 : USE cp_log_handling, ONLY: cp_get_default_logger,&
49 : cp_logger_get_default_io_unit,&
50 : cp_logger_get_default_unit_nr,&
51 : cp_logger_type,&
52 : cp_to_string
53 : USE cp_output_handling, ONLY: cp_p_file,&
54 : cp_print_key_should_output
55 : USE efield_utils, ONLY: efield_potential_lengh_gauge
56 : USE input_constants, ONLY: do_arnoldi,&
57 : do_bch,&
58 : do_em,&
59 : do_pade,&
60 : do_taylor
61 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
62 : section_vals_type
63 : USE iterate_matrix, ONLY: matrix_sqrt_Newton_Schulz
64 : USE kinds, ONLY: dp
65 : USE ls_matrix_exp, ONLY: cp_complex_dbcsr_gemm_3
66 : USE mathlib, ONLY: binomial
67 : USE parallel_gemm_api, ONLY: parallel_gemm
68 : USE particle_list_types, ONLY: particle_list_type
69 : USE pw_env_types, ONLY: pw_env_get,&
70 : pw_env_type
71 : USE pw_pool_types, ONLY: pw_pool_type
72 : USE pw_types, ONLY: pw_c1d_gs_type,&
73 : pw_r3d_rs_type
74 : USE qs_energy_init, ONLY: qs_energies_init
75 : USE qs_energy_types, ONLY: qs_energy_type
76 : USE qs_environment_types, ONLY: get_qs_env,&
77 : qs_environment_type
78 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
79 : USE qs_ks_types, ONLY: set_ks_env
80 : USE qs_loc_dipole, ONLY: loc_dipole
81 : USE qs_loc_states, ONLY: get_localization_info
82 : USE qs_loc_types, ONLY: qs_loc_env_create,&
83 : qs_loc_env_release,&
84 : qs_loc_env_type
85 : USE qs_loc_utils, ONLY: qs_loc_control_init,&
86 : qs_loc_init
87 : USE qs_mo_types, ONLY: get_mo_set,&
88 : mo_set_type
89 : USE rt_make_propagators, ONLY: propagate_arnoldi,&
90 : propagate_bch,&
91 : propagate_exp,&
92 : propagate_exp_density
93 : USE rt_propagation_output, ONLY: report_density_occupation,&
94 : rt_convergence,&
95 : rt_convergence_density
96 : USE rt_propagation_types, ONLY: get_rtp,&
97 : rt_prop_type
98 : USE rt_propagation_utils, ONLY: calc_S_derivs,&
99 : calc_update_rho,&
100 : calc_update_rho_sparse
101 : USE rt_propagation_velocity_gauge, ONLY: update_vector_potential,&
102 : velocity_gauge_ks_matrix
103 : #include "../base/base_uses.f90"
104 :
105 : IMPLICIT NONE
106 :
107 : PRIVATE
108 :
109 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_methods'
110 :
111 : PUBLIC :: propagation_step, &
112 : s_matrices_create, &
113 : calc_sinvH, &
114 : put_data_to_history, &
115 : rtp_localize
116 :
117 : CONTAINS
118 :
119 : ! **************************************************************************************************
120 : !> \brief performs a single propagation step a(t+Dt)=U(t+Dt,t)*a(0)
121 : !> and calculates the new exponential
122 : !> \param qs_env ...
123 : !> \param rtp ...
124 : !> \param rtp_control ...
125 : !> \author Florian Schiffmann (02.09)
126 : ! **************************************************************************************************
127 :
128 2316 : SUBROUTINE propagation_step(qs_env, rtp, rtp_control)
129 :
130 : TYPE(qs_environment_type), POINTER :: qs_env
131 : TYPE(rt_prop_type), POINTER :: rtp
132 : TYPE(rtp_control_type), POINTER :: rtp_control
133 :
134 : CHARACTER(len=*), PARAMETER :: routineN = 'propagation_step'
135 :
136 : INTEGER :: aspc_order, handle, i, im, re, unit_nr
137 : TYPE(cell_type), POINTER :: cell
138 2316 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: delta_mos, mos_new
139 : TYPE(cp_logger_type), POINTER :: logger
140 2316 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: delta_P, H_last_iter, ks_mix, ks_mix_im, &
141 2316 : matrix_ks, matrix_ks_im, matrix_s, &
142 2316 : rho_new
143 : TYPE(dft_control_type), POINTER :: dft_control
144 :
145 2316 : CALL timeset(routineN, handle)
146 :
147 2316 : logger => cp_get_default_logger()
148 2316 : IF (logger%para_env%is_source()) THEN
149 1158 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
150 : ELSE
151 : unit_nr = -1
152 : END IF
153 :
154 2316 : NULLIFY (cell, delta_P, rho_new, delta_mos, mos_new)
155 2316 : NULLIFY (ks_mix, ks_mix_im)
156 : ! get everything needed and set some values
157 2316 : CALL get_qs_env(qs_env, cell=cell, matrix_s=matrix_s, dft_control=dft_control)
158 :
159 2316 : IF (rtp%iter == 1) THEN
160 626 : CALL qs_energies_init(qs_env, .FALSE.)
161 : !the above recalculates matrix_s, but matrix not changed if ions are fixed
162 626 : IF (rtp_control%fixed_ions) CALL set_ks_env(qs_env%ks_env, s_mstruct_changed=.FALSE.)
163 :
164 : ! add additional terms to matrix_h and matrix_h_im in the case of applied electric field,
165 : ! either in the lengh or velocity gauge.
166 : ! should be called after qs_energies_init and before qs_ks_update_qs_env
167 626 : IF (dft_control%apply_efield_field) THEN
168 216 : IF (ANY(cell%perd(1:3) /= 0)) THEN
169 0 : CPABORT("Length gauge (efield) and periodicity are not compatible")
170 : END IF
171 54 : CALL efield_potential_lengh_gauge(qs_env)
172 572 : ELSE IF (rtp_control%velocity_gauge) THEN
173 32 : IF (dft_control%apply_vector_potential) THEN
174 32 : CALL update_vector_potential(qs_env, dft_control)
175 : END IF
176 32 : CALL velocity_gauge_ks_matrix(qs_env, subtract_nl_term=.FALSE.)
177 : END IF
178 :
179 626 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
180 626 : IF (.NOT. rtp_control%fixed_ions) THEN
181 274 : CALL s_matrices_create(matrix_s, rtp)
182 : END IF
183 626 : rtp%delta_iter = 100.0_dp
184 626 : rtp%mixing_factor = 1.0_dp
185 626 : rtp%mixing = .FALSE.
186 626 : aspc_order = rtp_control%aspc_order
187 626 : CALL aspc_extrapolate(rtp, matrix_s, aspc_order)
188 626 : IF (rtp%linear_scaling) THEN
189 182 : CALL calc_update_rho_sparse(qs_env)
190 : ELSE
191 444 : CALL calc_update_rho(qs_env)
192 : END IF
193 626 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE.)
194 : END IF
195 2316 : IF (.NOT. rtp_control%fixed_ions) THEN
196 1144 : CALL calc_S_derivs(qs_env)
197 : END IF
198 2316 : rtp%converged = .FALSE.
199 :
200 2316 : IF (rtp%linear_scaling) THEN
201 : ! keep temporary copy of the starting density matrix to check for convergence
202 806 : CALL get_rtp(rtp=rtp, rho_new=rho_new)
203 806 : NULLIFY (delta_P)
204 806 : CALL dbcsr_allocate_matrix_set(delta_P, SIZE(rho_new))
205 2910 : DO i = 1, SIZE(rho_new)
206 2104 : CALL dbcsr_init_p(delta_P(i)%matrix)
207 2104 : CALL dbcsr_create(delta_P(i)%matrix, template=rho_new(i)%matrix)
208 2910 : CALL dbcsr_copy(delta_P(i)%matrix, rho_new(i)%matrix)
209 : END DO
210 : ELSE
211 : ! keep temporary copy of the starting mos to check for convergence
212 1510 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
213 8390 : ALLOCATE (delta_mos(SIZE(mos_new)))
214 5370 : DO i = 1, SIZE(mos_new)
215 : CALL cp_fm_create(delta_mos(i), &
216 : matrix_struct=mos_new(i)%matrix_struct, &
217 3860 : name="delta_mos"//TRIM(ADJUSTL(cp_to_string(i))))
218 5370 : CALL cp_fm_to_fm(mos_new(i), delta_mos(i))
219 : END DO
220 : END IF
221 :
222 : CALL get_qs_env(qs_env, &
223 : matrix_ks=matrix_ks, &
224 2316 : matrix_ks_im=matrix_ks_im)
225 :
226 2316 : CALL get_rtp(rtp=rtp, H_last_iter=H_last_iter)
227 2316 : IF (rtp%mixing) THEN
228 96 : IF (unit_nr > 0) THEN
229 48 : WRITE (unit_nr, '(t3,a,2f16.8)') "Mixing the Hamiltonians to improve robustness, mixing factor: ", rtp%mixing_factor
230 : END IF
231 96 : CALL dbcsr_allocate_matrix_set(ks_mix, SIZE(matrix_ks))
232 96 : CALL dbcsr_allocate_matrix_set(ks_mix_im, SIZE(matrix_ks))
233 192 : DO i = 1, SIZE(matrix_ks)
234 96 : CALL dbcsr_init_p(ks_mix(i)%matrix)
235 96 : CALL dbcsr_create(ks_mix(i)%matrix, template=matrix_ks(1)%matrix)
236 96 : CALL dbcsr_init_p(ks_mix_im(i)%matrix)
237 192 : CALL dbcsr_create(ks_mix_im(i)%matrix, template=matrix_ks(1)%matrix, matrix_type=dbcsr_type_antisymmetric)
238 : END DO
239 192 : DO i = 1, SIZE(matrix_ks)
240 96 : re = 2*i - 1
241 96 : im = 2*i
242 96 : CALL dbcsr_add(ks_mix(i)%matrix, matrix_ks(i)%matrix, 0.0_dp, rtp%mixing_factor)
243 96 : CALL dbcsr_add(ks_mix(i)%matrix, H_last_iter(re)%matrix, 1.0_dp, 1.0_dp - rtp%mixing_factor)
244 192 : IF (rtp%propagate_complex_ks) THEN
245 0 : CALL dbcsr_add(ks_mix_im(i)%matrix, matrix_ks_im(i)%matrix, 0.0_dp, rtp%mixing_factor)
246 0 : CALL dbcsr_add(ks_mix_im(i)%matrix, H_last_iter(im)%matrix, 1.0_dp, 1.0_dp - rtp%mixing_factor)
247 : END IF
248 : END DO
249 96 : CALL calc_SinvH(rtp, ks_mix, ks_mix_im, rtp_control)
250 192 : DO i = 1, SIZE(matrix_ks)
251 96 : re = 2*i - 1
252 96 : im = 2*i
253 96 : CALL dbcsr_copy(H_last_iter(re)%matrix, ks_mix(i)%matrix)
254 192 : IF (rtp%propagate_complex_ks) THEN
255 0 : CALL dbcsr_copy(H_last_iter(im)%matrix, ks_mix_im(i)%matrix)
256 : END IF
257 : END DO
258 96 : CALL dbcsr_deallocate_matrix_set(ks_mix)
259 96 : CALL dbcsr_deallocate_matrix_set(ks_mix_im)
260 : ELSE
261 2220 : CALL calc_SinvH(rtp, matrix_ks, matrix_ks_im, rtp_control)
262 5106 : DO i = 1, SIZE(matrix_ks)
263 2886 : re = 2*i - 1
264 2886 : im = 2*i
265 2886 : CALL dbcsr_copy(H_last_iter(re)%matrix, matrix_ks(i)%matrix)
266 5106 : IF (rtp%propagate_complex_ks) THEN
267 438 : CALL dbcsr_copy(H_last_iter(im)%matrix, matrix_ks_im(i)%matrix)
268 : END IF
269 : END DO
270 : END IF
271 :
272 2316 : CALL compute_propagator_matrix(rtp, rtp_control%propagator)
273 :
274 4016 : SELECT CASE (rtp_control%mat_exp)
275 : CASE (do_pade, do_taylor)
276 1700 : IF (rtp%linear_scaling) THEN
277 622 : CALL propagate_exp_density(rtp, rtp_control)
278 622 : CALL calc_update_rho_sparse(qs_env)
279 : ELSE
280 1078 : CALL propagate_exp(rtp, rtp_control)
281 1078 : CALL calc_update_rho(qs_env)
282 : END IF
283 : CASE (do_arnoldi)
284 432 : CALL propagate_arnoldi(rtp, rtp_control)
285 432 : CALL calc_update_rho(qs_env)
286 : CASE (do_bch)
287 184 : CALL propagate_bch(rtp, rtp_control)
288 2500 : CALL calc_update_rho_sparse(qs_env)
289 : END SELECT
290 2316 : CALL step_finalize(qs_env, rtp_control, delta_mos, delta_P)
291 2316 : IF (rtp%linear_scaling) THEN
292 806 : CALL dbcsr_deallocate_matrix_set(delta_P)
293 : ELSE
294 1510 : CALL cp_fm_release(delta_mos)
295 : END IF
296 :
297 2316 : CALL timestop(handle)
298 :
299 2316 : END SUBROUTINE propagation_step
300 :
301 : ! **************************************************************************************************
302 : !> \brief Performs all the stuff to finish the step:
303 : !> convergence checks
304 : !> copying stuff into right place for the next step
305 : !> updating the history for extrapolation
306 : !> \param qs_env ...
307 : !> \param rtp_control ...
308 : !> \param delta_mos ...
309 : !> \param delta_P ...
310 : !> \author Florian Schiffmann (02.09)
311 : ! **************************************************************************************************
312 :
313 2316 : SUBROUTINE step_finalize(qs_env, rtp_control, delta_mos, delta_P)
314 : TYPE(qs_environment_type), POINTER :: qs_env
315 : TYPE(rtp_control_type), POINTER :: rtp_control
316 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: delta_mos
317 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: delta_P
318 :
319 : CHARACTER(len=*), PARAMETER :: routineN = 'step_finalize'
320 :
321 : INTEGER :: handle, i, ihist
322 2316 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new, mos_old
323 2316 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: exp_H_new, exp_H_old, matrix_ks, &
324 2316 : matrix_ks_im, rho_new, rho_old, s_mat
325 : TYPE(qs_energy_type), POINTER :: energy
326 : TYPE(rt_prop_type), POINTER :: rtp
327 :
328 2316 : CALL timeset(routineN, handle)
329 :
330 : CALL get_qs_env(qs_env=qs_env, rtp=rtp, matrix_s=s_mat, &
331 2316 : matrix_ks=matrix_ks, matrix_ks_im=matrix_ks_im, energy=energy)
332 2316 : CALL get_rtp(rtp=rtp, exp_H_old=exp_H_old, exp_H_new=exp_H_new)
333 :
334 2316 : IF (rtp_control%sc_check_start < rtp%iter) THEN
335 2316 : rtp%delta_iter_old = rtp%delta_iter
336 2316 : IF (rtp%linear_scaling) THEN
337 806 : CALL rt_convergence_density(rtp, delta_P, rtp%delta_iter)
338 : ELSE
339 1510 : CALL rt_convergence(rtp, s_mat(1)%matrix, delta_mos, rtp%delta_iter)
340 : END IF
341 2316 : rtp%converged = (rtp%delta_iter < rtp_control%eps_ener)
342 : !Apply mixing if scf loop is not converging
343 :
344 : !It would be better to redo the current step with mixixng,
345 : !but currently the decision is made to use mixing from the next step on
346 2316 : IF (rtp_control%sc_check_start < rtp%iter + 1) THEN
347 2316 : IF (rtp%delta_iter/rtp%delta_iter_old > 0.9) THEN
348 6 : rtp%mixing_factor = MAX(rtp%mixing_factor/2.0_dp, 0.125_dp)
349 6 : rtp%mixing = .TRUE.
350 : END IF
351 : END IF
352 : END IF
353 :
354 2316 : IF (rtp%converged) THEN
355 626 : IF (rtp%linear_scaling) THEN
356 182 : CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
357 : CALL purify_mcweeny_complex_nonorth(rho_new, s_mat, rtp%filter_eps, rtp%filter_eps_small, &
358 182 : rtp_control%mcweeny_max_iter, rtp_control%mcweeny_eps)
359 182 : IF (rtp_control%mcweeny_max_iter > 0) CALL calc_update_rho_sparse(qs_env)
360 182 : CALL report_density_occupation(rtp%filter_eps, rho_new)
361 686 : DO i = 1, SIZE(rho_new)
362 686 : CALL dbcsr_copy(rho_old(i)%matrix, rho_new(i)%matrix)
363 : END DO
364 : ELSE
365 444 : CALL get_rtp(rtp=rtp, mos_old=mos_old, mos_new=mos_new)
366 1544 : DO i = 1, SIZE(mos_new)
367 1544 : CALL cp_fm_to_fm(mos_new(i), mos_old(i))
368 : END DO
369 : END IF
370 626 : IF (rtp_control%propagator == do_em) CALL calc_SinvH(rtp, matrix_ks, matrix_ks_im, rtp_control)
371 2230 : DO i = 1, SIZE(exp_H_new)
372 2230 : CALL dbcsr_copy(exp_H_old(i)%matrix, exp_H_new(i)%matrix)
373 : END DO
374 626 : ihist = MOD(rtp%istep, rtp_control%aspc_order) + 1
375 626 : IF (rtp_control%fixed_ions) THEN
376 352 : CALL put_data_to_history(rtp, rho=rho_new, mos=mos_new, ihist=ihist)
377 : ELSE
378 274 : CALL put_data_to_history(rtp, rho=rho_new, mos=mos_new, s_mat=s_mat, ihist=ihist)
379 : END IF
380 : END IF
381 :
382 2316 : rtp%energy_new = energy%total
383 :
384 2316 : CALL timestop(handle)
385 :
386 2316 : END SUBROUTINE step_finalize
387 :
388 : ! **************************************************************************************************
389 : !> \brief computes the propagator matrix for EM/ETRS, RTP/EMD
390 : !> \param rtp ...
391 : !> \param propagator ...
392 : !> \author Florian Schiffmann (02.09)
393 : ! **************************************************************************************************
394 :
395 4632 : SUBROUTINE compute_propagator_matrix(rtp, propagator)
396 : TYPE(rt_prop_type), POINTER :: rtp
397 : INTEGER :: propagator
398 :
399 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_propagator_matrix'
400 :
401 : INTEGER :: handle, i
402 : REAL(Kind=dp) :: dt, prefac
403 2316 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: exp_H_new, exp_H_old, propagator_matrix
404 :
405 2316 : CALL timeset(routineN, handle)
406 : CALL get_rtp(rtp=rtp, exp_H_new=exp_H_new, exp_H_old=exp_H_old, &
407 2316 : propagator_matrix=propagator_matrix, dt=dt)
408 :
409 2316 : prefac = -0.5_dp*dt
410 :
411 8280 : DO i = 1, SIZE(exp_H_new)
412 5964 : CALL dbcsr_add(propagator_matrix(i)%matrix, exp_H_new(i)%matrix, 0.0_dp, prefac)
413 8280 : IF (propagator == do_em) THEN
414 64 : CALL dbcsr_add(propagator_matrix(i)%matrix, exp_H_old(i)%matrix, 1.0_dp, prefac)
415 : END IF
416 : END DO
417 :
418 2316 : CALL timestop(handle)
419 :
420 2316 : END SUBROUTINE compute_propagator_matrix
421 :
422 : ! **************************************************************************************************
423 : !> \brief computes S_inv*H, if needed Sinv*B and S_inv*H_imag and store these quantities to the
424 : !> \brief exp_H for the real and imag part (for RTP and EMD)
425 : !> \param rtp ...
426 : !> \param matrix_ks ...
427 : !> \param matrix_ks_im ...
428 : !> \param rtp_control ...
429 : !> \author Florian Schiffmann (02.09)
430 : ! **************************************************************************************************
431 :
432 2532 : SUBROUTINE calc_SinvH(rtp, matrix_ks, matrix_ks_im, rtp_control)
433 : TYPE(rt_prop_type), POINTER :: rtp
434 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_im
435 : TYPE(rtp_control_type), POINTER :: rtp_control
436 :
437 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_SinvH'
438 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
439 :
440 : INTEGER :: handle, im, ispin, re
441 2532 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: exp_H, SinvB, SinvH, SinvH_imag
442 : TYPE(dbcsr_type) :: matrix_ks_nosym
443 : TYPE(dbcsr_type), POINTER :: B_mat, S_inv
444 :
445 2532 : CALL timeset(routineN, handle)
446 2532 : CALL get_rtp(rtp=rtp, S_inv=S_inv, exp_H_new=exp_H)
447 5800 : DO ispin = 1, SIZE(matrix_ks)
448 3268 : re = ispin*2 - 1
449 3268 : im = ispin*2
450 3268 : CALL dbcsr_set(exp_H(re)%matrix, zero)
451 5800 : CALL dbcsr_set(exp_H(im)%matrix, zero)
452 : END DO
453 2532 : CALL dbcsr_create(matrix_ks_nosym, template=matrix_ks(1)%matrix, matrix_type="N")
454 :
455 : ! Real part of S_inv x H -> imag part of exp_H
456 5800 : DO ispin = 1, SIZE(matrix_ks)
457 3268 : re = ispin*2 - 1
458 3268 : im = ispin*2
459 3268 : CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, matrix_ks_nosym)
460 : CALL dbcsr_multiply("N", "N", one, S_inv, matrix_ks_nosym, zero, exp_H(im)%matrix, &
461 3268 : filter_eps=rtp%filter_eps)
462 5800 : IF (.NOT. rtp_control%fixed_ions) THEN
463 1478 : CALL get_rtp(rtp=rtp, SinvH=SinvH)
464 1478 : CALL dbcsr_copy(SinvH(ispin)%matrix, exp_H(im)%matrix)
465 : END IF
466 : END DO
467 :
468 : ! Imag part of S_inv x H -> real part of exp_H
469 2532 : IF (rtp%propagate_complex_ks) THEN
470 940 : DO ispin = 1, SIZE(matrix_ks)
471 486 : re = ispin*2 - 1
472 486 : im = ispin*2
473 486 : CALL dbcsr_set(matrix_ks_nosym, 0.0_dp)
474 486 : CALL dbcsr_desymmetrize(matrix_ks_im(ispin)%matrix, matrix_ks_nosym)
475 : ! - SinvH_imag is added to exp_H(re)%matrix
476 : CALL dbcsr_multiply("N", "N", -one, S_inv, matrix_ks_nosym, zero, exp_H(re)%matrix, &
477 486 : filter_eps=rtp%filter_eps)
478 940 : IF (.NOT. rtp_control%fixed_ions) THEN
479 286 : CALL get_rtp(rtp=rtp, SinvH_imag=SinvH_imag)
480 : ! -SinvH_imag is saved
481 286 : CALL dbcsr_copy(SinvH_imag(ispin)%matrix, exp_H(re)%matrix)
482 : END IF
483 : END DO
484 : END IF
485 : ! EMD case: the real part of exp_H should be updated with B
486 2532 : IF (.NOT. rtp_control%fixed_ions) THEN
487 1226 : CALL get_rtp(rtp=rtp, B_mat=B_mat, SinvB=SinvB)
488 1226 : CALL dbcsr_set(matrix_ks_nosym, 0.0_dp)
489 1226 : CALL dbcsr_multiply("N", "N", one, S_inv, B_mat, zero, matrix_ks_nosym, filter_eps=rtp%filter_eps)
490 2704 : DO ispin = 1, SIZE(matrix_ks)
491 1478 : re = ispin*2 - 1
492 1478 : im = ispin*2
493 : ! + SinvB is added to exp_H(re)%matrix
494 1478 : CALL dbcsr_add(exp_H(re)%matrix, matrix_ks_nosym, 1.0_dp, 1.0_dp)
495 : ! + SinvB is saved
496 2704 : CALL dbcsr_copy(SinvB(ispin)%matrix, matrix_ks_nosym)
497 : END DO
498 : END IF
499 : ! Otherwise no real part for exp_H
500 :
501 2532 : CALL dbcsr_release(matrix_ks_nosym)
502 2532 : CALL timestop(handle)
503 :
504 2532 : END SUBROUTINE calc_SinvH
505 :
506 : ! **************************************************************************************************
507 : !> \brief calculates the needed overlap-like matrices
508 : !> depending on the way the exponential is calculated, only S^-1 is needed
509 : !> \param s_mat ...
510 : !> \param rtp ...
511 : !> \author Florian Schiffmann (02.09)
512 : ! **************************************************************************************************
513 :
514 478 : SUBROUTINE s_matrices_create(s_mat, rtp)
515 :
516 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: s_mat
517 : TYPE(rt_prop_type), POINTER :: rtp
518 :
519 : CHARACTER(len=*), PARAMETER :: routineN = 's_matrices_create'
520 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
521 :
522 : INTEGER :: handle
523 : TYPE(dbcsr_type), POINTER :: S_half, S_inv, S_minus_half
524 :
525 478 : CALL timeset(routineN, handle)
526 :
527 478 : CALL get_rtp(rtp=rtp, S_inv=S_inv)
528 :
529 478 : IF (rtp%linear_scaling) THEN
530 136 : CALL get_rtp(rtp=rtp, S_half=S_half, S_minus_half=S_minus_half)
531 : CALL matrix_sqrt_Newton_Schulz(S_half, S_minus_half, s_mat(1)%matrix, rtp%filter_eps, &
532 136 : rtp%newton_schulz_order, rtp%lanzcos_threshold, rtp%lanzcos_max_iter)
533 : CALL dbcsr_multiply("N", "N", one, S_minus_half, S_minus_half, zero, S_inv, &
534 136 : filter_eps=rtp%filter_eps)
535 : ELSE
536 342 : CALL dbcsr_copy(S_inv, s_mat(1)%matrix)
537 : CALL cp_dbcsr_cholesky_decompose(S_inv, para_env=rtp%ao_ao_fmstruct%para_env, &
538 342 : blacs_env=rtp%ao_ao_fmstruct%context)
539 : CALL cp_dbcsr_cholesky_invert(S_inv, para_env=rtp%ao_ao_fmstruct%para_env, &
540 342 : blacs_env=rtp%ao_ao_fmstruct%context, uplo_to_full=.TRUE.)
541 : END IF
542 :
543 478 : CALL timestop(handle)
544 478 : END SUBROUTINE s_matrices_create
545 :
546 : ! **************************************************************************************************
547 : !> \brief Calculates the frobenius norm of a complex matrix represented by two real matrices
548 : !> \param frob_norm ...
549 : !> \param mat_re ...
550 : !> \param mat_im ...
551 : !> \author Samuel Andermatt (04.14)
552 : ! **************************************************************************************************
553 :
554 568 : SUBROUTINE complex_frobenius_norm(frob_norm, mat_re, mat_im)
555 :
556 : REAL(KIND=dp), INTENT(out) :: frob_norm
557 : TYPE(dbcsr_type), POINTER :: mat_re, mat_im
558 :
559 : CHARACTER(len=*), PARAMETER :: routineN = 'complex_frobenius_norm'
560 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
561 :
562 : INTEGER :: col_atom, handle, row_atom
563 : LOGICAL :: found
564 284 : REAL(dp), DIMENSION(:, :), POINTER :: block_values, block_values2
565 : TYPE(dbcsr_iterator_type) :: iter
566 : TYPE(dbcsr_type), POINTER :: tmp
567 :
568 284 : CALL timeset(routineN, handle)
569 :
570 : NULLIFY (tmp)
571 284 : ALLOCATE (tmp)
572 284 : CALL dbcsr_create(tmp, template=mat_re)
573 : !make sure the tmp has the same sparsity pattern as the real and the complex part combined
574 284 : CALL dbcsr_add(tmp, mat_re, zero, one)
575 284 : CALL dbcsr_add(tmp, mat_im, zero, one)
576 284 : CALL dbcsr_set(tmp, zero)
577 : !calculate the hadamard product
578 284 : CALL dbcsr_iterator_start(iter, tmp)
579 1804 : DO WHILE (dbcsr_iterator_blocks_left(iter))
580 1520 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
581 1520 : CALL dbcsr_get_block_p(mat_re, row_atom, col_atom, block_values2, found=found)
582 1520 : IF (found) THEN
583 264788 : block_values = block_values2*block_values2
584 : END IF
585 1520 : CALL dbcsr_get_block_p(mat_im, row_atom, col_atom, block_values2, found=found)
586 1520 : IF (found) THEN
587 250308 : block_values = block_values + block_values2*block_values2
588 : END IF
589 133438 : block_values = SQRT(block_values)
590 : END DO
591 284 : CALL dbcsr_iterator_stop(iter)
592 284 : frob_norm = dbcsr_frobenius_norm(tmp)
593 :
594 284 : CALL dbcsr_deallocate_matrix(tmp)
595 :
596 284 : CALL timestop(handle)
597 :
598 284 : END SUBROUTINE complex_frobenius_norm
599 :
600 : ! **************************************************************************************************
601 : !> \brief Does McWeeny for complex matrices in the non-orthogonal basis
602 : !> \param P ...
603 : !> \param s_mat ...
604 : !> \param eps ...
605 : !> \param eps_small ...
606 : !> \param max_iter ...
607 : !> \param threshold ...
608 : !> \author Samuel Andermatt (04.14)
609 : ! **************************************************************************************************
610 :
611 182 : SUBROUTINE purify_mcweeny_complex_nonorth(P, s_mat, eps, eps_small, max_iter, threshold)
612 :
613 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: P, s_mat
614 : REAL(KIND=dp), INTENT(in) :: eps, eps_small
615 : INTEGER, INTENT(in) :: max_iter
616 : REAL(KIND=dp), INTENT(in) :: threshold
617 :
618 : CHARACTER(len=*), PARAMETER :: routineN = 'purify_mcweeny_complex_nonorth'
619 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
620 :
621 : INTEGER :: handle, i, im, imax, ispin, re, unit_nr
622 : REAL(KIND=dp) :: frob_norm
623 : TYPE(cp_logger_type), POINTER :: logger
624 182 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: PS, PSP, tmp
625 :
626 182 : CALL timeset(routineN, handle)
627 :
628 182 : logger => cp_get_default_logger()
629 182 : IF (logger%para_env%is_source()) THEN
630 91 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
631 : ELSE
632 : unit_nr = -1
633 : END IF
634 :
635 182 : NULLIFY (tmp, PS, PSP)
636 182 : CALL dbcsr_allocate_matrix_set(tmp, SIZE(P))
637 182 : CALL dbcsr_allocate_matrix_set(PSP, SIZE(P))
638 182 : CALL dbcsr_allocate_matrix_set(PS, SIZE(P))
639 686 : DO i = 1, SIZE(P)
640 504 : CALL dbcsr_init_p(PS(i)%matrix)
641 504 : CALL dbcsr_create(PS(i)%matrix, template=P(1)%matrix)
642 504 : CALL dbcsr_init_p(PSP(i)%matrix)
643 504 : CALL dbcsr_create(PSP(i)%matrix, template=P(1)%matrix)
644 504 : CALL dbcsr_init_p(tmp(i)%matrix)
645 686 : CALL dbcsr_create(tmp(i)%matrix, template=P(1)%matrix)
646 : END DO
647 182 : IF (SIZE(P) == 2) THEN
648 112 : CALL dbcsr_scale(P(1)%matrix, one/2)
649 112 : CALL dbcsr_scale(P(2)%matrix, one/2)
650 : END IF
651 434 : DO ispin = 1, SIZE(P)/2
652 252 : re = 2*ispin - 1
653 252 : im = 2*ispin
654 252 : imax = MAX(max_iter, 1) !if max_iter is 0 then only the deviation from idempotency needs to be calculated
655 516 : DO i = 1, imax
656 : CALL dbcsr_multiply("N", "N", one, P(re)%matrix, s_mat(1)%matrix, &
657 284 : zero, PS(re)%matrix, filter_eps=eps_small)
658 : CALL dbcsr_multiply("N", "N", one, P(im)%matrix, s_mat(1)%matrix, &
659 284 : zero, PS(im)%matrix, filter_eps=eps_small)
660 : CALL cp_complex_dbcsr_gemm_3("N", "N", one, PS(re)%matrix, PS(im)%matrix, &
661 : P(re)%matrix, P(im)%matrix, zero, PSP(re)%matrix, PSP(im)%matrix, &
662 284 : filter_eps=eps_small)
663 284 : CALL dbcsr_copy(tmp(re)%matrix, PSP(re)%matrix)
664 284 : CALL dbcsr_copy(tmp(im)%matrix, PSP(im)%matrix)
665 284 : CALL dbcsr_add(tmp(re)%matrix, P(re)%matrix, 1.0_dp, -1.0_dp)
666 284 : CALL dbcsr_add(tmp(im)%matrix, P(im)%matrix, 1.0_dp, -1.0_dp)
667 284 : CALL complex_frobenius_norm(frob_norm, tmp(re)%matrix, tmp(im)%matrix)
668 284 : IF (unit_nr > 0) WRITE (unit_nr, '(t3,a,2f16.8)') "Deviation from idempotency: ", frob_norm
669 800 : IF (frob_norm > threshold .AND. max_iter > 0) THEN
670 264 : CALL dbcsr_copy(P(re)%matrix, PSP(re)%matrix)
671 264 : CALL dbcsr_copy(P(im)%matrix, PSP(im)%matrix)
672 : CALL cp_complex_dbcsr_gemm_3("N", "N", -2.0_dp, PS(re)%matrix, PS(im)%matrix, &
673 : PSP(re)%matrix, PSP(im)%matrix, 3.0_dp, P(re)%matrix, P(im)%matrix, &
674 264 : filter_eps=eps_small)
675 264 : CALL dbcsr_filter(P(re)%matrix, eps)
676 264 : CALL dbcsr_filter(P(im)%matrix, eps)
677 : !make sure P is exactly hermitian
678 264 : CALL dbcsr_transposed(tmp(re)%matrix, P(re)%matrix)
679 264 : CALL dbcsr_add(P(re)%matrix, tmp(re)%matrix, one/2, one/2)
680 264 : CALL dbcsr_transposed(tmp(im)%matrix, P(im)%matrix)
681 264 : CALL dbcsr_add(P(im)%matrix, tmp(im)%matrix, one/2, -one/2)
682 : ELSE
683 : EXIT
684 : END IF
685 : END DO
686 : !make sure P is hermitian
687 252 : CALL dbcsr_transposed(tmp(re)%matrix, P(re)%matrix)
688 252 : CALL dbcsr_add(P(re)%matrix, tmp(re)%matrix, one/2, one/2)
689 252 : CALL dbcsr_transposed(tmp(im)%matrix, P(im)%matrix)
690 434 : CALL dbcsr_add(P(im)%matrix, tmp(im)%matrix, one/2, -one/2)
691 : END DO
692 182 : IF (SIZE(P) == 2) THEN
693 112 : CALL dbcsr_scale(P(1)%matrix, one*2)
694 112 : CALL dbcsr_scale(P(2)%matrix, one*2)
695 : END IF
696 182 : CALL dbcsr_deallocate_matrix_set(tmp)
697 182 : CALL dbcsr_deallocate_matrix_set(PS)
698 182 : CALL dbcsr_deallocate_matrix_set(PSP)
699 :
700 182 : CALL timestop(handle)
701 :
702 182 : END SUBROUTINE purify_mcweeny_complex_nonorth
703 :
704 : ! **************************************************************************************************
705 : !> \brief ...
706 : !> \param rtp ...
707 : !> \param matrix_s ...
708 : !> \param aspc_order ...
709 : ! **************************************************************************************************
710 626 : SUBROUTINE aspc_extrapolate(rtp, matrix_s, aspc_order)
711 : TYPE(rt_prop_type), POINTER :: rtp
712 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
713 : INTEGER, INTENT(in) :: aspc_order
714 :
715 : CHARACTER(len=*), PARAMETER :: routineN = 'aspc_extrapolate'
716 : COMPLEX(KIND=dp), PARAMETER :: cone = (1.0_dp, 0.0_dp), &
717 : czero = (0.0_dp, 0.0_dp)
718 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
719 :
720 : INTEGER :: handle, i, iaspc, icol_local, ihist, &
721 : imat, k, kdbl, n, naspc, ncol_local, &
722 : nmat
723 : REAL(KIND=dp) :: alpha
724 : TYPE(cp_cfm_type) :: cfm_tmp, cfm_tmp1, csc
725 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct, matrix_struct_new
726 : TYPE(cp_fm_type) :: fm_tmp, fm_tmp1, fm_tmp2
727 626 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
728 626 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: mo_hist
729 626 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new, s_hist
730 626 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_hist
731 :
732 626 : NULLIFY (rho_hist)
733 626 : CALL timeset(routineN, handle)
734 626 : CALL cite_reference(Kolafa2004)
735 626 : CALL cite_reference(Kuhne2007)
736 :
737 626 : IF (rtp%linear_scaling) THEN
738 182 : CALL get_rtp(rtp=rtp, rho_new=rho_new)
739 : ELSE
740 444 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
741 : END IF
742 :
743 626 : naspc = MIN(rtp%istep, aspc_order)
744 626 : IF (rtp%linear_scaling) THEN
745 182 : nmat = SIZE(rho_new)
746 182 : rho_hist => rtp%history%rho_history
747 686 : DO imat = 1, nmat
748 1454 : DO iaspc = 1, naspc
749 : alpha = (-1.0_dp)**(iaspc + 1)*REAL(iaspc, KIND=dp)* &
750 768 : binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1)
751 768 : ihist = MOD(rtp%istep - iaspc, aspc_order) + 1
752 1272 : IF (iaspc == 1) THEN
753 504 : CALL dbcsr_add(rho_new(imat)%matrix, rho_hist(imat, ihist)%matrix, zero, alpha)
754 : ELSE
755 264 : CALL dbcsr_add(rho_new(imat)%matrix, rho_hist(imat, ihist)%matrix, one, alpha)
756 : END IF
757 : END DO
758 : END DO
759 : ELSE
760 444 : mo_hist => rtp%history%mo_history
761 444 : nmat = SIZE(mos_new)
762 1544 : DO imat = 1, nmat
763 3968 : DO iaspc = 1, naspc
764 : alpha = (-1.0_dp)**(iaspc + 1)*REAL(iaspc, KIND=dp)* &
765 2424 : binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1)
766 2424 : ihist = MOD(rtp%istep - iaspc, aspc_order) + 1
767 3524 : IF (iaspc == 1) THEN
768 1100 : CALL cp_fm_scale_and_add(zero, mos_new(imat), alpha, mo_hist(imat, ihist))
769 : ELSE
770 1324 : CALL cp_fm_scale_and_add(one, mos_new(imat), alpha, mo_hist(imat, ihist))
771 : END IF
772 : END DO
773 : END DO
774 :
775 444 : mo_hist => rtp%history%mo_history
776 444 : s_hist => rtp%history%s_history
777 994 : DO i = 1, SIZE(mos_new)/2
778 550 : NULLIFY (matrix_struct, matrix_struct_new)
779 :
780 : CALL cp_fm_struct_double(matrix_struct, &
781 : mos_new(2*i)%matrix_struct, &
782 : mos_new(2*i)%matrix_struct%context, &
783 550 : .TRUE., .FALSE.)
784 :
785 550 : CALL cp_fm_create(fm_tmp, matrix_struct)
786 550 : CALL cp_fm_create(fm_tmp1, matrix_struct)
787 550 : CALL cp_fm_create(fm_tmp2, mos_new(2*i)%matrix_struct)
788 550 : CALL cp_cfm_create(cfm_tmp, mos_new(2*i)%matrix_struct)
789 550 : CALL cp_cfm_create(cfm_tmp1, mos_new(2*i)%matrix_struct)
790 :
791 550 : CALL cp_fm_get_info(fm_tmp, ncol_global=kdbl)
792 :
793 : CALL cp_fm_get_info(mos_new(2*i), &
794 : nrow_global=n, &
795 : ncol_global=k, &
796 550 : ncol_local=ncol_local)
797 :
798 : CALL cp_fm_struct_create(matrix_struct_new, &
799 : template_fmstruct=mos_new(2*i)%matrix_struct, &
800 : nrow_global=k, &
801 550 : ncol_global=k)
802 550 : CALL cp_cfm_create(csc, matrix_struct_new)
803 :
804 550 : CALL cp_fm_struct_release(matrix_struct_new)
805 550 : CALL cp_fm_struct_release(matrix_struct)
806 :
807 : ! first the most recent
808 :
809 : ! reorthogonalize vectors
810 2612 : DO icol_local = 1, ncol_local
811 22188 : fm_tmp%local_data(:, icol_local) = mos_new(2*i - 1)%local_data(:, icol_local)
812 22738 : fm_tmp%local_data(:, icol_local + ncol_local) = mos_new(2*i)%local_data(:, icol_local)
813 : END DO
814 :
815 550 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, fm_tmp, fm_tmp1, kdbl)
816 :
817 2612 : DO icol_local = 1, ncol_local
818 : cfm_tmp%local_data(:, icol_local) = CMPLX(fm_tmp1%local_data(:, icol_local), &
819 22188 : fm_tmp1%local_data(:, icol_local + ncol_local), dp)
820 : cfm_tmp1%local_data(:, icol_local) = CMPLX(mos_new(2*i - 1)%local_data(:, icol_local), &
821 22738 : mos_new(2*i)%local_data(:, icol_local), dp)
822 : END DO
823 550 : CALL parallel_gemm('C', 'N', k, k, n, cone, cfm_tmp1, cfm_tmp, czero, csc)
824 550 : CALL cp_cfm_cholesky_decompose(csc)
825 550 : CALL cp_cfm_triangular_multiply(csc, cfm_tmp1, n_cols=k, side='R', invert_tr=.TRUE.)
826 2612 : DO icol_local = 1, ncol_local
827 22188 : mos_new(2*i - 1)%local_data(:, icol_local) = REAL(cfm_tmp1%local_data(:, icol_local), dp)
828 22738 : mos_new(2*i)%local_data(:, icol_local) = AIMAG(cfm_tmp1%local_data(:, icol_local))
829 : END DO
830 :
831 : ! deallocate work matrices
832 550 : CALL cp_cfm_release(csc)
833 550 : CALL cp_fm_release(fm_tmp)
834 550 : CALL cp_fm_release(fm_tmp1)
835 550 : CALL cp_fm_release(fm_tmp2)
836 550 : CALL cp_cfm_release(cfm_tmp)
837 2094 : CALL cp_cfm_release(cfm_tmp1)
838 : END DO
839 :
840 : END IF
841 :
842 626 : CALL timestop(handle)
843 :
844 626 : END SUBROUTINE aspc_extrapolate
845 :
846 : ! **************************************************************************************************
847 : !> \brief ...
848 : !> \param rtp ...
849 : !> \param mos ...
850 : !> \param rho ...
851 : !> \param s_mat ...
852 : !> \param ihist ...
853 : ! **************************************************************************************************
854 830 : SUBROUTINE put_data_to_history(rtp, mos, rho, s_mat, ihist)
855 : TYPE(rt_prop_type), POINTER :: rtp
856 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos
857 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho
858 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
859 : POINTER :: s_mat
860 : INTEGER :: ihist
861 :
862 : INTEGER :: i
863 :
864 830 : IF (rtp%linear_scaling) THEN
865 1032 : DO i = 1, SIZE(rho)
866 1032 : CALL dbcsr_copy(rtp%history%rho_history(i, ihist)%matrix, rho(i)%matrix)
867 : END DO
868 : ELSE
869 1950 : DO i = 1, SIZE(mos)
870 1950 : CALL cp_fm_to_fm(mos(i), rtp%history%mo_history(i, ihist))
871 : END DO
872 558 : IF (PRESENT(s_mat)) THEN
873 342 : IF (ASSOCIATED(rtp%history%s_history(ihist)%matrix)) THEN ! the sparsity might be different
874 : ! (future struct:check)
875 124 : CALL dbcsr_deallocate_matrix(rtp%history%s_history(ihist)%matrix)
876 : END IF
877 342 : ALLOCATE (rtp%history%s_history(ihist)%matrix)
878 342 : CALL dbcsr_copy(rtp%history%s_history(ihist)%matrix, s_mat(1)%matrix)
879 : END IF
880 : END IF
881 :
882 830 : END SUBROUTINE put_data_to_history
883 :
884 : ! **************************************************************************************************
885 : !> \brief Computes Maximally localised Wannier functions and print properties according to
886 : !> FORCE_EVAL%DFT%LOCALIZE, adapted from qs_scf_post_gpw::scf_post_calculation_gpw
887 : !> \param qs_env QuickStep environment
888 : !> \param rtp Real time propagation environment
889 : !> \par History 03/2020 created [LS]
890 : !> \author Lukas Schreder
891 : ! **************************************************************************************************
892 352 : SUBROUTINE rtp_localize(qs_env, rtp)
893 :
894 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
895 : TYPE(rt_prop_type), INTENT(IN), POINTER :: rtp
896 :
897 : CHARACTER(len=*), PARAMETER :: routineN = 'rtp_localize'
898 :
899 : INTEGER :: handle, ispin, output_unit
900 352 : INTEGER, DIMENSION(:, :, :), POINTER :: marked_states
901 : LOGICAL :: do_homo, do_mo_cubes, do_wannier_cubes, &
902 : p_loc
903 352 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
904 352 : TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: occupied_evals
905 352 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mo_localized, occupied_orbs
906 352 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_coeff
907 : TYPE(cp_logger_type), POINTER :: logger
908 : TYPE(dft_control_type), POINTER :: dft_control
909 352 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
910 : TYPE(particle_list_type), POINTER :: particles
911 : TYPE(pw_c1d_gs_type) :: wf_g
912 : TYPE(pw_env_type), POINTER :: pw_env
913 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
914 : TYPE(pw_r3d_rs_type) :: wf_r
915 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
916 : TYPE(section_vals_type), POINTER :: dft_section, input, loc_print_section, &
917 : loc_section, print_key
918 :
919 352 : CALL timeset(routineN, handle)
920 :
921 : ! Localization of propagated orbitals requires MO coefficients
922 352 : IF (rtp%linear_scaling) THEN
923 136 : CALL timestop(handle)
924 272 : RETURN
925 : END IF
926 :
927 216 : CALL cite_reference(Schreder2021)
928 :
929 216 : NULLIFY (auxbas_pw_pool, dft_control, dft_section, input, loc_print_section, &
930 216 : loc_section, logger, marked_states, mo_coeff, mo_eigenvalues, mos, &
931 216 : occupied_evals, particles, print_key, pw_env, qs_loc_env)
932 :
933 216 : logger => cp_get_default_logger()
934 216 : output_unit = cp_logger_get_default_io_unit(logger)
935 :
936 216 : IF (output_unit > 0) THEN
937 108 : WRITE (unit=output_unit, fmt="(A)") "LOCALIZE| Localizing propagated orbitals"
938 : END IF
939 :
940 : ! get section properties
941 216 : CALL get_qs_env(qs_env, dft_control=dft_control, input=input, pw_env=pw_env)
942 : ! get propagated MO coeffs
943 216 : CALL get_rtp(rtp, mos_new=mo_coeff)
944 216 : dft_section => section_vals_get_subs_vals(input, "DFT")
945 216 : loc_section => section_vals_get_subs_vals(dft_section, "LOCALIZE")
946 216 : loc_print_section => section_vals_get_subs_vals(loc_section, "PRINT")
947 :
948 : ! what properties to print out
949 216 : print_key => section_vals_get_subs_vals(loc_print_section, "MOLECULAR_DIPOLES")
950 216 : p_loc = BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
951 :
952 216 : print_key => section_vals_get_subs_vals(loc_print_section, "TOTAL_DIPOLE")
953 : p_loc = p_loc &
954 216 : .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
955 216 : print_key => section_vals_get_subs_vals(loc_print_section, "WANNIER_CENTERS")
956 : p_loc = p_loc &
957 216 : .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
958 216 : print_key => section_vals_get_subs_vals(loc_print_section, "WANNIER_SPREADS")
959 : p_loc = p_loc &
960 216 : .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
961 216 : print_key => section_vals_get_subs_vals(loc_print_section, "WANNIER_CUBES")
962 : p_loc = p_loc &
963 216 : .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
964 216 : print_key => section_vals_get_subs_vals(loc_print_section, "MOLECULAR_STATES")
965 : p_loc = p_loc &
966 216 : .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
967 216 : print_key => section_vals_get_subs_vals(loc_print_section, "MOLECULAR_MOMENTS")
968 : p_loc = p_loc &
969 216 : .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
970 216 : print_key => section_vals_get_subs_vals(loc_print_section, "LOCALIZED_MOMENTS")
971 : p_loc = p_loc &
972 216 : .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
973 216 : print_key => section_vals_get_subs_vals(loc_print_section, "WANNIER_STATES")
974 : p_loc = p_loc &
975 216 : .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
976 :
977 : do_wannier_cubes = BTEST(cp_print_key_should_output(logger%iter_info, loc_print_section, &
978 216 : "WANNIER_CUBES"), cp_p_file)
979 :
980 216 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
981 216 : CALL auxbas_pw_pool%create_pw(wf_r)
982 216 : CALL auxbas_pw_pool%create_pw(wf_g)
983 :
984 216 : IF (p_loc) THEN
985 122 : ALLOCATE (occupied_evals(dft_control%nspins))
986 166 : ALLOCATE (occupied_orbs(SIZE(mo_coeff)))
987 140 : ALLOCATE (mo_localized(SIZE(mo_coeff)))
988 26 : CALL get_qs_env(qs_env, mos=mos)
989 26 : CALL get_rtp(rtp, mos_new=mo_coeff)
990 114 : DO ispin = 1, SIZE(mo_coeff)
991 88 : occupied_orbs(ispin) = mo_coeff(ispin)
992 88 : CALL cp_fm_create(mo_localized(ispin), mo_coeff(ispin)%matrix_struct)
993 114 : CALL cp_fm_to_fm(mo_coeff(ispin), mo_localized(ispin))
994 : END DO
995 :
996 70 : DO ispin = 1, dft_control%nspins
997 44 : CALL get_mo_set(mos(ispin), eigenvalues=mo_eigenvalues)
998 70 : occupied_evals(ispin)%array => mo_eigenvalues
999 : END DO
1000 :
1001 26 : do_homo = .TRUE.
1002 182 : ALLOCATE (qs_loc_env)
1003 26 : CALL qs_loc_env_create(qs_loc_env)
1004 26 : CALL qs_loc_control_init(qs_loc_env, loc_section, do_homo=do_homo)
1005 26 : CALL qs_loc_init(qs_env, qs_loc_env, loc_section, mo_localized, do_homo, do_mo_cubes)
1006 : CALL get_localization_info(qs_env, qs_loc_env, loc_section, mo_localized, wf_r, wf_g, &
1007 26 : particles, occupied_orbs, occupied_evals, marked_states)
1008 26 : CALL loc_dipole(input, dft_control, qs_loc_env, logger, qs_env)
1009 :
1010 114 : DO ispin = 1, SIZE(mo_localized)
1011 114 : CALL cp_fm_release(mo_localized(ispin))
1012 : END DO
1013 26 : DEALLOCATE (mo_localized)
1014 26 : DEALLOCATE (occupied_orbs)
1015 26 : DEALLOCATE (occupied_evals)
1016 26 : CALL qs_loc_env_release(qs_loc_env)
1017 26 : DEALLOCATE (qs_loc_env)
1018 52 : IF (ASSOCIATED(marked_states)) THEN
1019 18 : DEALLOCATE (marked_states)
1020 : END IF
1021 : END IF
1022 :
1023 216 : CALL auxbas_pw_pool%give_back_pw(wf_r)
1024 216 : CALL auxbas_pw_pool%give_back_pw(wf_g)
1025 :
1026 216 : CALL timestop(handle)
1027 :
1028 840 : END SUBROUTINE rtp_localize
1029 :
1030 : END MODULE rt_propagation_methods
|
