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