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