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 needed for EMD
10 : !> \author Florian Schiffmann (02.09)
11 : ! **************************************************************************************************
12 :
13 : MODULE rt_propagation_utils
14 : USE atomic_kind_types, ONLY: atomic_kind_type
15 : USE cell_types, ONLY: cell_type
16 : USE cp_blacs_env, ONLY: cp_blacs_env_type
17 : USE cp_control_types, ONLY: dft_control_type,&
18 : rtp_control_type
19 : USE cp_dbcsr_api, ONLY: &
20 : dbcsr_add, dbcsr_binary_read, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, &
21 : dbcsr_desymmetrize, dbcsr_distribution_type, dbcsr_filter, dbcsr_get_info, &
22 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
23 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_scale, dbcsr_set, dbcsr_type
24 : USE cp_dbcsr_contrib, ONLY: dbcsr_checksum
25 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t,&
26 : dbcsr_deallocate_matrix_set
27 : USE cp_files, ONLY: close_file,&
28 : file_exists,&
29 : open_file
30 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
31 : USE cp_fm_types, ONLY: cp_fm_create,&
32 : cp_fm_get_info,&
33 : cp_fm_release,&
34 : cp_fm_set_all,&
35 : cp_fm_to_fm,&
36 : cp_fm_type
37 : USE cp_log_handling, ONLY: cp_get_default_logger,&
38 : cp_logger_get_default_io_unit,&
39 : cp_logger_get_default_unit_nr,&
40 : cp_logger_type
41 : USE cp_output_handling, ONLY: cp_p_file,&
42 : cp_print_key_finished_output,&
43 : cp_print_key_generate_filename,&
44 : cp_print_key_should_output,&
45 : cp_print_key_unit_nr
46 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
47 : USE efield_utils, ONLY: make_field
48 : USE input_constants, ONLY: use_restart_wfn,&
49 : use_rt_restart
50 : USE input_section_types, ONLY: section_get_ival,&
51 : section_get_ivals,&
52 : section_get_lval,&
53 : section_vals_get,&
54 : section_vals_get_subs_vals,&
55 : section_vals_type,&
56 : section_vals_val_get
57 : USE kinds, ONLY: default_path_length,&
58 : default_string_length,&
59 : dp
60 : USE mathconstants, ONLY: zero
61 : USE memory_utilities, ONLY: reallocate
62 : USE message_passing, ONLY: mp_para_env_type
63 : USE orbital_pointers, ONLY: ncoset
64 : USE particle_list_types, ONLY: particle_list_type
65 : USE particle_types, ONLY: particle_type
66 : USE physcon, ONLY: femtoseconds
67 : USE pw_env_types, ONLY: pw_env_get,&
68 : pw_env_type
69 : USE pw_methods, ONLY: pw_multiply,&
70 : pw_zero
71 : USE pw_pool_types, ONLY: pw_pool_p_type,&
72 : pw_pool_type
73 : USE pw_types, ONLY: pw_c1d_gs_type,&
74 : pw_r3d_rs_type
75 : USE qs_collocate_density, ONLY: calculate_wavefunction
76 : USE qs_density_matrices, ONLY: calculate_density_matrix
77 : USE qs_dftb_matrices, ONLY: build_dftb_overlap
78 : USE qs_environment_types, ONLY: get_qs_env,&
79 : qs_environment_type
80 : USE qs_kind_types, ONLY: qs_kind_type
81 : USE qs_ks_types, ONLY: qs_ks_did_change,&
82 : qs_ks_env_type
83 : USE qs_mo_io, ONLY: read_mo_set_from_restart,&
84 : read_rt_mos_from_restart,&
85 : write_mo_set_to_output_unit
86 : USE qs_mo_types, ONLY: allocate_mo_set,&
87 : deallocate_mo_set,&
88 : get_mo_set,&
89 : init_mo_set,&
90 : mo_set_type
91 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
92 : USE qs_overlap, ONLY: build_overlap_matrix
93 : USE qs_rho_methods, ONLY: qs_rho_update_rho
94 : USE qs_rho_types, ONLY: qs_rho_get,&
95 : qs_rho_set,&
96 : qs_rho_type
97 : USE qs_scf_wfn_mix, ONLY: wfn_mix
98 : USE qs_subsys_types, ONLY: qs_subsys_get,&
99 : qs_subsys_type
100 : USE rt_propagation_types, ONLY: get_rtp,&
101 : rt_prop_type
102 : #include "../base/base_uses.f90"
103 :
104 : IMPLICIT NONE
105 : PRIVATE
106 :
107 : PUBLIC :: get_restart_wfn, &
108 : calc_S_derivs, &
109 : calc_update_rho, &
110 : calc_update_rho_sparse, &
111 : calculate_P_imaginary, &
112 : write_rtp_mos_to_output_unit, &
113 : write_rtp_mo_cubes, &
114 : warn_section_unused, &
115 : read_moments, &
116 : recalculate_fields
117 :
118 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_utils'
119 :
120 : CONTAINS
121 :
122 : ! **************************************************************************************************
123 : !> \brief Calculates dS/dR respectily the velocity weighted derivatves
124 : !> only needed for ehrenfest MD.
125 : !>
126 : !> \param qs_env the qs environment
127 : !> \par History
128 : !> 02.2009 created [Manuel Guidon]
129 : !> 02.2014 switched to dbcsr matrices [Samuel Andermatt]
130 : !> \author Florian Schiffmann
131 : ! **************************************************************************************************
132 1218 : SUBROUTINE calc_S_derivs(qs_env)
133 : TYPE(qs_environment_type), POINTER :: qs_env
134 :
135 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_S_derivs'
136 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
137 :
138 : INTEGER :: col_atom, handle, i, j, m, maxder, n, &
139 : nder, row_atom
140 : INTEGER, DIMENSION(6, 2) :: c_map_mat
141 : LOGICAL :: return_s_derivatives
142 1218 : REAL(dp), DIMENSION(:, :), POINTER :: block_values
143 : TYPE(dbcsr_iterator_type) :: iter
144 1218 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: C_mat, S_der, s_derivs
145 : TYPE(dbcsr_type), POINTER :: B_mat, tmp_mat, tmp_mat2
146 : TYPE(dft_control_type), POINTER :: dft_control
147 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
148 1218 : POINTER :: sab_orb
149 1218 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
150 : TYPE(qs_ks_env_type), POINTER :: ks_env
151 : TYPE(rt_prop_type), POINTER :: rtp
152 :
153 1218 : CALL timeset(routineN, handle)
154 :
155 1218 : return_s_derivatives = .TRUE.
156 :
157 1218 : NULLIFY (particle_set)
158 1218 : NULLIFY (rtp)
159 1218 : NULLIFY (s_derivs)
160 1218 : NULLIFY (dft_control)
161 1218 : NULLIFY (ks_env)
162 :
163 : CALL get_qs_env(qs_env=qs_env, &
164 : rtp=rtp, &
165 : particle_set=particle_set, &
166 : sab_orb=sab_orb, &
167 : dft_control=dft_control, &
168 1218 : ks_env=ks_env)
169 :
170 1218 : CALL get_rtp(rtp=rtp, B_mat=B_mat, C_mat=C_mat, S_der=S_der)
171 :
172 1218 : nder = 2
173 1218 : maxder = ncoset(nder)
174 :
175 : NULLIFY (tmp_mat)
176 1218 : ALLOCATE (tmp_mat)
177 1218 : CALL dbcsr_create(tmp_mat, template=S_der(1)%matrix, matrix_type="N")
178 :
179 1218 : IF (rtp%iter < 2) THEN
180 : ! calculate the overlap derivative matrices
181 348 : IF (dft_control%qs_control%dftb) THEN
182 84 : CALL build_dftb_overlap(qs_env, nder, s_derivs)
183 : ELSE
184 : CALL build_overlap_matrix(ks_env, nderivative=nder, matrix_s=s_derivs, &
185 264 : basis_type_a="ORB", basis_type_b="ORB", sab_nl=sab_orb)
186 : END IF
187 :
188 : NULLIFY (tmp_mat2)
189 348 : ALLOCATE (tmp_mat2)
190 348 : CALL dbcsr_create(tmp_mat2, template=S_der(1)%matrix, matrix_type="S")
191 3480 : DO m = 1, 9
192 3132 : CALL dbcsr_copy(tmp_mat2, s_derivs(m + 1)%matrix)
193 3132 : CALL dbcsr_desymmetrize(tmp_mat2, S_der(m)%matrix)
194 3132 : CALL dbcsr_scale(S_der(m)%matrix, -one)
195 3132 : CALL dbcsr_filter(S_der(m)%matrix, rtp%filter_eps)
196 : !The diagonal should be zero
197 3132 : CALL dbcsr_iterator_start(iter, S_der(m)%matrix)
198 15335 : DO WHILE (dbcsr_iterator_blocks_left(iter))
199 12203 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
200 194291 : IF (row_atom == col_atom) block_values = 0.0_dp
201 : END DO
202 6612 : CALL dbcsr_iterator_stop(iter)
203 : END DO
204 348 : CALL dbcsr_deallocate_matrix_set(s_derivs)
205 348 : CALL dbcsr_deallocate_matrix(tmp_mat2)
206 : END IF
207 :
208 : !calculate scalar product v(Rb)*<alpha|d/dRb beta> (B_mat), and store the first derivatives
209 :
210 1218 : CALL dbcsr_set(B_mat, zero)
211 4872 : DO m = 1, 3
212 3654 : CALL dbcsr_copy(tmp_mat, S_der(m)%matrix)
213 3654 : CALL dbcsr_iterator_start(iter, tmp_mat)
214 18766 : DO WHILE (dbcsr_iterator_blocks_left(iter))
215 15112 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
216 222322 : IF (row_atom == col_atom) block_values = 0.0_dp
217 466536 : block_values = block_values*particle_set(col_atom)%v(m)
218 : END DO
219 3654 : CALL dbcsr_iterator_stop(iter)
220 8526 : CALL dbcsr_add(B_mat, tmp_mat, one, one)
221 : END DO
222 1218 : CALL dbcsr_filter(B_mat, rtp%filter_eps)
223 : !calculate C matrix: v(Rb)*<d/dRa alpha| d/dRb beta>
224 :
225 1218 : c_map_mat = 0
226 1218 : n = 0
227 4872 : DO j = 1, 3
228 12180 : DO m = j, 3
229 7308 : n = n + 1
230 7308 : c_map_mat(n, 1) = j
231 7308 : IF (m == j) CYCLE
232 10962 : c_map_mat(n, 2) = m
233 : END DO
234 : END DO
235 :
236 4872 : DO i = 1, 3
237 4872 : CALL dbcsr_set(C_mat(i)%matrix, zero)
238 : END DO
239 8526 : DO m = 1, 6
240 7308 : CALL dbcsr_copy(tmp_mat, S_der(m + 3)%matrix)
241 23142 : DO j = 1, 2
242 14616 : IF (c_map_mat(m, j) == 0) CYCLE
243 21924 : CALL dbcsr_add(C_mat(c_map_mat(m, j))%matrix, tmp_mat, one, one)
244 : END DO
245 : END DO
246 :
247 4872 : DO m = 1, 3
248 3654 : CALL dbcsr_iterator_start(iter, C_mat(m)%matrix)
249 17716 : DO WHILE (dbcsr_iterator_blocks_left(iter))
250 14062 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
251 459585 : block_values = block_values*particle_set(row_atom)%v(m)
252 : END DO
253 3654 : CALL dbcsr_iterator_stop(iter)
254 8526 : CALL dbcsr_filter(C_mat(m)%matrix, rtp%filter_eps)
255 : END DO
256 :
257 1218 : CALL dbcsr_deallocate_matrix(tmp_mat)
258 1218 : CALL timestop(handle)
259 1218 : END SUBROUTINE calc_S_derivs
260 :
261 : ! **************************************************************************************************
262 : !> \brief reads the restart file. At the moment only SCF (means only real)
263 : !> \param qs_env ...
264 : !> \author Florian Schiffmann (02.09)
265 : ! **************************************************************************************************
266 :
267 36 : SUBROUTINE get_restart_wfn(qs_env)
268 : TYPE(qs_environment_type), POINTER :: qs_env
269 :
270 : CHARACTER(LEN=default_path_length) :: file_name, project_name
271 : INTEGER :: i, id_nr, im, ispin, ncol, nspin, &
272 : output_unit, re, unit_nr
273 : REAL(KIND=dp) :: alpha, cs_pos
274 36 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
275 : TYPE(cp_fm_type) :: mos_occ
276 36 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_old
277 : TYPE(cp_logger_type), POINTER :: logger
278 : TYPE(dbcsr_distribution_type) :: dist
279 36 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv, rho_new, rho_old
280 : TYPE(dft_control_type), POINTER :: dft_control
281 36 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
282 : TYPE(mp_para_env_type), POINTER :: para_env
283 36 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
284 36 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
285 : TYPE(qs_rho_type), POINTER :: rho_struct
286 : TYPE(rt_prop_type), POINTER :: rtp
287 : TYPE(section_vals_type), POINTER :: dft_section, input
288 :
289 36 : NULLIFY (atomic_kind_set, qs_kind_set, mo_array, particle_set, rho_struct, para_env)
290 :
291 : CALL get_qs_env(qs_env, &
292 : qs_kind_set=qs_kind_set, &
293 : atomic_kind_set=atomic_kind_set, &
294 : particle_set=particle_set, &
295 : mos=mo_array, &
296 : input=input, &
297 : rtp=rtp, &
298 : dft_control=dft_control, &
299 : rho=rho_struct, &
300 36 : para_env=para_env)
301 36 : logger => cp_get_default_logger()
302 36 : output_unit = cp_logger_get_default_io_unit(logger)
303 :
304 36 : IF (logger%para_env%is_source()) THEN
305 18 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
306 : ELSE
307 : unit_nr = -1
308 : END IF
309 :
310 36 : id_nr = 0
311 36 : nspin = SIZE(mo_array)
312 36 : CALL qs_rho_get(rho_struct, rho_ao=p_rmpv)
313 36 : dft_section => section_vals_get_subs_vals(input, "DFT")
314 62 : SELECT CASE (dft_control%rtp_control%initial_wfn)
315 : CASE (use_restart_wfn)
316 : CALL read_mo_set_from_restart(mo_array, qs_kind_set, particle_set, para_env, &
317 : id_nr=id_nr, multiplicity=dft_control%multiplicity, &
318 26 : dft_section=dft_section)
319 26 : CALL set_uniform_occupation_mo_array(mo_array, nspin)
320 :
321 26 : IF (dft_control%rtp_control%apply_wfn_mix_init_restart) THEN
322 : CALL wfn_mix(mo_array, particle_set, dft_section, qs_kind_set, para_env, output_unit, &
323 4 : for_rtp=.TRUE.)
324 : END IF
325 :
326 70 : DO ispin = 1, nspin
327 70 : CALL calculate_density_matrix(mo_array(ispin), p_rmpv(ispin)%matrix)
328 : END DO
329 26 : IF (rtp%linear_scaling) THEN
330 14 : CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
331 34 : DO ispin = 1, nspin
332 20 : re = 2*ispin - 1
333 20 : im = 2*ispin
334 20 : CALL cp_fm_get_info(mo_array(ispin)%mo_coeff, ncol_global=ncol)
335 : CALL cp_fm_create(mos_occ, &
336 : matrix_struct=mo_array(ispin)%mo_coeff%matrix_struct, &
337 20 : name="mos_occ")
338 20 : CALL cp_fm_to_fm(mo_array(ispin)%mo_coeff, mos_occ)
339 20 : IF (mo_array(ispin)%uniform_occupation) THEN
340 16 : alpha = 3.0_dp - REAL(nspin, dp)
341 78 : CALL cp_fm_column_scale(mos_occ, mo_array(ispin)%occupation_numbers/alpha)
342 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
343 : matrix_v=mos_occ, &
344 : ncol=ncol, &
345 16 : alpha=alpha, keep_sparsity=.FALSE.)
346 : ELSE
347 4 : alpha = 1.0_dp
348 88 : CALL cp_fm_column_scale(mos_occ, mo_array(ispin)%occupation_numbers/alpha)
349 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
350 : matrix_v=mo_array(ispin)%mo_coeff, &
351 : matrix_g=mos_occ, &
352 : ncol=ncol, &
353 4 : alpha=alpha, keep_sparsity=.FALSE.)
354 : END IF
355 20 : CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
356 20 : CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
357 54 : CALL cp_fm_release(mos_occ)
358 : END DO
359 14 : CALL calc_update_rho_sparse(qs_env)
360 : ELSE
361 12 : CALL get_rtp(rtp=rtp, mos_old=mos_old)
362 36 : DO i = 1, SIZE(qs_env%mos)
363 24 : CALL cp_fm_to_fm(mo_array(i)%mo_coeff, mos_old(2*i - 1))
364 36 : CALL cp_fm_set_all(mos_old(2*i), zero, zero)
365 : END DO
366 : END IF
367 : CASE (use_rt_restart)
368 36 : IF (rtp%linear_scaling) THEN
369 2 : CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
370 2 : project_name = logger%iter_info%project_name
371 4 : DO ispin = 1, nspin
372 2 : re = 2*ispin - 1
373 2 : im = 2*ispin
374 2 : WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_RE", ispin, "_RESTART.dm"
375 2 : CALL dbcsr_get_info(rho_old(re)%matrix, distribution=dist)
376 2 : CALL dbcsr_binary_read(file_name, distribution=dist, matrix_new=rho_old(re)%matrix)
377 2 : cs_pos = dbcsr_checksum(rho_old(re)%matrix, pos=.TRUE.)
378 2 : IF (unit_nr > 0) THEN
379 1 : WRITE (unit_nr, '(T2,A,E20.8)') "Read restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
380 : END IF
381 2 : WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_IM", ispin, "_RESTART.dm"
382 2 : CALL dbcsr_get_info(rho_old(im)%matrix, distribution=dist)
383 2 : CALL dbcsr_binary_read(file_name, distribution=dist, matrix_new=rho_old(im)%matrix)
384 2 : cs_pos = dbcsr_checksum(rho_old(im)%matrix, pos=.TRUE.)
385 8 : IF (unit_nr > 0) THEN
386 1 : WRITE (unit_nr, '(T2,A,E20.8)') "Read restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
387 : END IF
388 : END DO
389 6 : DO i = 1, SIZE(rho_new)
390 6 : CALL dbcsr_copy(rho_new(i)%matrix, rho_old(i)%matrix)
391 : END DO
392 2 : CALL calc_update_rho_sparse(qs_env)
393 : ELSE
394 8 : CALL get_rtp(rtp=rtp, mos_old=mos_old)
395 : CALL read_rt_mos_from_restart(mo_array, mos_old, qs_kind_set, particle_set, para_env, &
396 8 : id_nr, dft_control%multiplicity, dft_section)
397 8 : CALL set_uniform_occupation_mo_array(mo_array, nspin)
398 16 : DO ispin = 1, nspin
399 : CALL calculate_density_matrix(mo_array(ispin), &
400 16 : p_rmpv(ispin)%matrix)
401 : END DO
402 : END IF
403 : END SELECT
404 :
405 36 : END SUBROUTINE get_restart_wfn
406 :
407 : ! **************************************************************************************************
408 : !> \brief Set mo_array(ispin)%uniform_occupation after a restart
409 : !> \param mo_array ...
410 : !> \param nspin ...
411 : !> \author Guillaume Le Breton (03.23)
412 : ! **************************************************************************************************
413 :
414 34 : SUBROUTINE set_uniform_occupation_mo_array(mo_array, nspin)
415 :
416 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
417 : INTEGER :: nspin
418 :
419 : INTEGER :: ispin, mo
420 : LOGICAL :: is_uniform
421 :
422 86 : DO ispin = 1, nspin
423 52 : is_uniform = .TRUE.
424 274 : DO mo = 1, mo_array(ispin)%nmo
425 : IF (mo_array(ispin)%occupation_numbers(mo) /= 0.0 .AND. &
426 222 : mo_array(ispin)%occupation_numbers(mo) /= 1.0 .AND. &
427 52 : mo_array(ispin)%occupation_numbers(mo) /= 2.0) THEN
428 48 : is_uniform = .FALSE.
429 : END IF
430 : END DO
431 86 : mo_array(ispin)%uniform_occupation = is_uniform
432 : END DO
433 :
434 34 : END SUBROUTINE set_uniform_occupation_mo_array
435 :
436 : ! **************************************************************************************************
437 : !> \brief calculates the density from the complex MOs and passes the density to
438 : !> qs_env.
439 : !> \param qs_env ...
440 : !> \author Florian Schiffmann (02.09)
441 : ! **************************************************************************************************
442 :
443 2068 : SUBROUTINE calc_update_rho(qs_env)
444 :
445 : TYPE(qs_environment_type), POINTER :: qs_env
446 :
447 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_update_rho'
448 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
449 :
450 : INTEGER :: handle, i, im, ncol, re
451 : REAL(KIND=dp) :: alpha
452 : TYPE(cp_fm_type) :: mos_occ
453 2068 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
454 2068 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao, rho_ao_im
455 2068 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
456 : TYPE(qs_ks_env_type), POINTER :: ks_env
457 : TYPE(qs_rho_type), POINTER :: rho
458 : TYPE(rt_prop_type), POINTER :: rtp
459 :
460 2068 : CALL timeset(routineN, handle)
461 :
462 2068 : NULLIFY (rho, ks_env, mos_new, rtp)
463 : CALL get_qs_env(qs_env, &
464 : ks_env=ks_env, &
465 : rho=rho, &
466 : rtp=rtp, &
467 2068 : mos=mos)
468 2068 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
469 2068 : CALL qs_rho_get(rho_struct=rho, rho_ao=rho_ao)
470 4694 : DO i = 1, SIZE(mos_new)/2
471 2626 : re = 2*i - 1; im = 2*i
472 2626 : CALL dbcsr_set(rho_ao(i)%matrix, zero)
473 2626 : CALL cp_fm_get_info(mos_new(re), ncol_global=ncol)
474 : CALL cp_fm_create(mos_occ, &
475 : matrix_struct=mos(i)%mo_coeff%matrix_struct, &
476 2626 : name="mos_occ")
477 2626 : CALL cp_fm_to_fm(mos_new(re), mos_occ)
478 2626 : IF (mos(i)%uniform_occupation) THEN
479 2526 : alpha = 3*one - REAL(SIZE(mos_new)/2, dp)
480 11926 : CALL cp_fm_column_scale(mos_occ, mos(i)%occupation_numbers/alpha)
481 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(i)%matrix, &
482 : matrix_v=mos_occ, &
483 : ncol=ncol, &
484 2526 : alpha=alpha)
485 : ELSE
486 100 : alpha = 1.0_dp
487 660 : CALL cp_fm_column_scale(mos_occ, mos(i)%occupation_numbers/alpha)
488 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(i)%matrix, &
489 : matrix_v=mos_new(re), &
490 : matrix_g=mos_occ, &
491 : ncol=ncol, &
492 100 : alpha=alpha)
493 : END IF
494 :
495 : ! It is actually complex conjugate but i*i=-1 therefore it must be added
496 2626 : CALL cp_fm_to_fm(mos_new(im), mos_occ)
497 2626 : IF (mos(i)%uniform_occupation) THEN
498 2526 : alpha = 3*one - REAL(SIZE(mos_new)/2, dp)
499 11926 : CALL cp_fm_column_scale(mos_occ, mos(i)%occupation_numbers/alpha)
500 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(i)%matrix, &
501 : matrix_v=mos_occ, &
502 : ncol=ncol, &
503 2526 : alpha=alpha)
504 : ELSE
505 100 : alpha = 1.0_dp
506 660 : CALL cp_fm_column_scale(mos_occ, mos(i)%occupation_numbers/alpha)
507 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(i)%matrix, &
508 : matrix_v=mos_new(im), &
509 : matrix_g=mos_occ, &
510 : ncol=ncol, &
511 100 : alpha=alpha)
512 : END IF
513 7320 : CALL cp_fm_release(mos_occ)
514 : END DO
515 2068 : CALL qs_rho_update_rho(rho, qs_env)
516 :
517 2068 : IF (rtp%track_imag_density) THEN
518 1408 : CALL qs_rho_get(rho_struct=rho, rho_ao_im=rho_ao_im)
519 1408 : CALL calculate_P_imaginary(qs_env, rtp, rho_ao_im, keep_sparsity=.TRUE.)
520 1408 : CALL qs_rho_set(rho, rho_ao_im=rho_ao_im)
521 : END IF
522 :
523 2068 : CALL qs_ks_did_change(ks_env, rho_changed=.TRUE.)
524 :
525 2068 : CALL timestop(handle)
526 :
527 2068 : END SUBROUTINE calc_update_rho
528 :
529 : ! **************************************************************************************************
530 : !> \brief Copies the density matrix back into the qs_env%rho%rho_ao
531 : !> \param qs_env ...
532 : !> \author Samuel Andermatt (3.14)
533 : ! **************************************************************************************************
534 :
535 1254 : SUBROUTINE calc_update_rho_sparse(qs_env)
536 :
537 : TYPE(qs_environment_type), POINTER :: qs_env
538 :
539 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_update_rho_sparse'
540 : REAL(KIND=dp), PARAMETER :: zero = 0.0_dp
541 :
542 : INTEGER :: handle, ispin
543 1254 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao, rho_ao_im, rho_new
544 : TYPE(dft_control_type), POINTER :: dft_control
545 : TYPE(qs_ks_env_type), POINTER :: ks_env
546 : TYPE(qs_rho_type), POINTER :: rho
547 : TYPE(rt_prop_type), POINTER :: rtp
548 : TYPE(rtp_control_type), POINTER :: rtp_control
549 :
550 1254 : NULLIFY (rho, ks_env, rtp, dft_control)
551 1254 : CALL timeset(routineN, handle)
552 : CALL get_qs_env(qs_env, &
553 : ks_env=ks_env, &
554 : rho=rho, &
555 : rtp=rtp, &
556 1254 : dft_control=dft_control)
557 1254 : rtp_control => dft_control%rtp_control
558 1254 : CALL get_rtp(rtp=rtp, rho_new=rho_new)
559 1254 : CALL qs_rho_get(rho_struct=rho, rho_ao=rho_ao)
560 1254 : IF (rtp%track_imag_density) CALL qs_rho_get(rho_struct=rho, rho_ao_im=rho_ao_im)
561 2926 : DO ispin = 1, SIZE(rho_ao)
562 1672 : CALL dbcsr_set(rho_ao(ispin)%matrix, zero)
563 1672 : CALL dbcsr_copy(rho_ao(ispin)%matrix, rho_new(ispin*2 - 1)%matrix, keep_sparsity=.TRUE.)
564 2926 : IF (rtp%track_imag_density) THEN
565 482 : CALL dbcsr_copy(rho_ao_im(ispin)%matrix, rho_new(ispin*2)%matrix, keep_sparsity=.TRUE.)
566 : END IF
567 : END DO
568 :
569 1254 : CALL qs_rho_update_rho(rho, qs_env)
570 1254 : CALL qs_ks_did_change(ks_env, rho_changed=.TRUE.)
571 :
572 1254 : CALL timestop(handle)
573 :
574 1254 : END SUBROUTINE calc_update_rho_sparse
575 :
576 : ! **************************************************************************************************
577 : !> \brief ...
578 : !> \param qs_env ...
579 : !> \param rtp ...
580 : !> \param matrix_p_im ...
581 : !> \param keep_sparsity ...
582 : ! **************************************************************************************************
583 1408 : SUBROUTINE calculate_P_imaginary(qs_env, rtp, matrix_p_im, keep_sparsity)
584 : TYPE(qs_environment_type), POINTER :: qs_env
585 : TYPE(rt_prop_type), POINTER :: rtp
586 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p_im
587 : LOGICAL, OPTIONAL :: keep_sparsity
588 :
589 : INTEGER :: i, im, ncol, re
590 : LOGICAL :: my_keep_sparsity
591 : REAL(KIND=dp) :: alpha
592 1408 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new, mos_occ
593 1408 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
594 :
595 1408 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
596 :
597 1408 : my_keep_sparsity = .FALSE.
598 1408 : IF (PRESENT(keep_sparsity)) my_keep_sparsity = keep_sparsity
599 1408 : CALL get_qs_env(qs_env, mos=mos)
600 7700 : ALLOCATE (mos_occ(SIZE(mos)*2))
601 :
602 3146 : DO i = 1, SIZE(mos_new)/2
603 1738 : re = 2*i - 1; im = 2*i
604 1738 : alpha = 3.0_dp - REAL(SIZE(matrix_p_im), dp)
605 : CALL cp_fm_create(mos_occ(re), &
606 : matrix_struct=mos(i)%mo_coeff%matrix_struct, &
607 1738 : name="mos_occ")
608 : CALL cp_fm_create(mos_occ(im), &
609 : matrix_struct=mos(i)%mo_coeff%matrix_struct, &
610 1738 : name="mos_occ")
611 1738 : CALL dbcsr_set(matrix_p_im(i)%matrix, 0.0_dp)
612 1738 : CALL cp_fm_get_info(mos_new(re), ncol_global=ncol)
613 1738 : CALL cp_fm_to_fm(mos_new(re), mos_occ(re))
614 8768 : CALL cp_fm_column_scale(mos_occ(re), mos(i)%occupation_numbers/alpha)
615 1738 : CALL cp_fm_to_fm(mos_new(im), mos_occ(im))
616 8768 : CALL cp_fm_column_scale(mos_occ(im), mos(i)%occupation_numbers/alpha)
617 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=matrix_p_im(i)%matrix, &
618 : matrix_v=mos_occ(im), &
619 : matrix_g=mos_occ(re), &
620 : ncol=ncol, &
621 : keep_sparsity=my_keep_sparsity, &
622 : alpha=2.0_dp*alpha, &
623 4884 : symmetry_mode=-1)
624 : END DO
625 1408 : CALL cp_fm_release(mos_occ)
626 :
627 1408 : END SUBROUTINE calculate_P_imaginary
628 :
629 : ! **************************************************************************************************
630 : !> \brief ...
631 : !> \param qs_env ...
632 : !> \param rtp ...
633 : ! **************************************************************************************************
634 648 : SUBROUTINE write_rtp_mos_to_output_unit(qs_env, rtp)
635 : TYPE(qs_environment_type), POINTER :: qs_env
636 : TYPE(rt_prop_type), POINTER :: rtp
637 :
638 : CHARACTER(len=*), PARAMETER :: routineN = 'write_rtp_mos_to_output_unit'
639 :
640 : CHARACTER(LEN=10) :: spin
641 : CHARACTER(LEN=2*default_string_length) :: name
642 : INTEGER :: handle, i, ispin, nao, nelectron, nmo, &
643 : nspins
644 : LOGICAL :: print_eigvecs, print_mo_info
645 : REAL(KIND=dp) :: flexible_electron_count, maxocc, n_el_f
646 324 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
647 324 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
648 : TYPE(cp_logger_type), POINTER :: logger
649 : TYPE(mo_set_type) :: mo_set_rtp
650 324 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
651 324 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
652 324 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
653 : TYPE(section_vals_type), POINTER :: dft_section, input
654 :
655 324 : NULLIFY (atomic_kind_set, particle_set, qs_kind_set, input, mos, dft_section)
656 :
657 324 : CALL timeset(routineN, handle)
658 :
659 : CALL get_qs_env(qs_env, &
660 : atomic_kind_set=atomic_kind_set, &
661 : qs_kind_set=qs_kind_set, &
662 : particle_set=particle_set, &
663 : input=input, &
664 324 : mos=mos)
665 : ! Quick return, if no printout of MO information is requested
666 324 : dft_section => section_vals_get_subs_vals(input, "DFT")
667 324 : CALL section_vals_val_get(dft_section, "PRINT%MO%EIGENVECTORS", l_val=print_eigvecs)
668 :
669 324 : NULLIFY (logger)
670 324 : logger => cp_get_default_logger()
671 : print_mo_info = (cp_print_key_should_output(logger%iter_info, &
672 : dft_section, "PRINT%MO") /= 0) .OR. &
673 324 : (qs_env%sim_step == 1)
674 :
675 106 : IF ((.NOT. print_mo_info) .OR. (.NOT. print_eigvecs)) THEN
676 320 : CALL timestop(handle)
677 320 : RETURN
678 : END IF
679 :
680 4 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
681 :
682 4 : nspins = SIZE(mos_new)/2
683 :
684 8 : DO ispin = 1, nspins
685 : ! initiate mo_set
686 : CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo, nelectron=nelectron, &
687 4 : n_el_f=n_el_f, maxocc=maxocc, flexible_electron_count=flexible_electron_count)
688 :
689 : CALL allocate_mo_set(mo_set_rtp, &
690 : nao=nao, &
691 : nmo=nmo, &
692 : nelectron=nelectron, &
693 : n_el_f=n_el_f, &
694 : maxocc=maxocc, &
695 4 : flexible_electron_count=flexible_electron_count)
696 :
697 4 : WRITE (name, FMT="(A,I6)") "RTP MO SET, SPIN ", ispin
698 4 : CALL init_mo_set(mo_set_rtp, fm_ref=mos_new(2*ispin - 1), name=name)
699 :
700 4 : IF (nspins > 1) THEN
701 0 : IF (ispin == 1) THEN
702 0 : spin = "ALPHA SPIN"
703 : ELSE
704 0 : spin = "BETA SPIN"
705 : END IF
706 : ELSE
707 4 : spin = ""
708 : END IF
709 :
710 8 : mo_set_rtp%occupation_numbers = mos(ispin)%occupation_numbers
711 :
712 : !loop for real (odd) and imaginary (even) parts
713 12 : DO i = 1, 0, -1
714 8 : CALL cp_fm_to_fm(mos_new(2*ispin - i), mo_set_rtp%mo_coeff)
715 : CALL write_mo_set_to_output_unit(mo_set_rtp, qs_kind_set, particle_set, &
716 : dft_section, 4, 0, rtp=.TRUE., spin=TRIM(spin), &
717 12 : cpart=MOD(i, 2), sim_step=qs_env%sim_step)
718 : END DO
719 :
720 12 : CALL deallocate_mo_set(mo_set_rtp)
721 : END DO
722 :
723 4 : CALL timestop(handle)
724 :
725 324 : END SUBROUTINE write_rtp_mos_to_output_unit
726 :
727 : ! **************************************************************************************************
728 : !> \brief Write the time dependent amplitude of the MOs in real grid.
729 : !> Very close to qs_scf_post_gpw/qs_scf_post_occ_cubes subroutine.
730 : !> \param qs_env ...
731 : !> \param rtp ...
732 : !> \author Guillaume Le Breton (11.22)
733 : ! **************************************************************************************************
734 324 : SUBROUTINE write_rtp_mo_cubes(qs_env, rtp)
735 : TYPE(qs_environment_type), POINTER :: qs_env
736 : TYPE(rt_prop_type), POINTER :: rtp
737 :
738 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_rtp_mo_cubes'
739 :
740 : CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title
741 : INTEGER :: handle, homo, i, ir, ispin, ivector, &
742 : n_rep, nhomo, nlist, nspins, &
743 : rt_time_step, unit_nr
744 324 : INTEGER, DIMENSION(:), POINTER :: list, list_index
745 : LOGICAL :: append_cube, do_kpoints, mpi_io
746 324 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
747 : TYPE(cell_type), POINTER :: cell
748 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
749 324 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
750 : TYPE(cp_fm_type), POINTER :: mo_coeff
751 : TYPE(cp_logger_type), POINTER :: logger
752 : TYPE(dft_control_type), POINTER :: dft_control
753 324 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
754 : TYPE(mp_para_env_type), POINTER :: para_env
755 : TYPE(particle_list_type), POINTER :: particles
756 324 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
757 : TYPE(pw_c1d_gs_type) :: wf_g
758 : TYPE(pw_env_type), POINTER :: pw_env
759 324 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
760 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
761 : TYPE(pw_r3d_rs_type) :: density_r, wf_r
762 324 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
763 : TYPE(qs_subsys_type), POINTER :: subsys
764 : TYPE(section_vals_type), POINTER :: dft_section, input
765 :
766 324 : CALL timeset(routineN, handle)
767 :
768 324 : NULLIFY (logger, auxbas_pw_pool, pw_pools, pw_env)
769 :
770 : ! Get all the info from qs:
771 : CALL get_qs_env(qs_env, do_kpoints=do_kpoints, &
772 324 : input=input)
773 :
774 : ! Kill the run in the case of K points
775 324 : IF (do_kpoints) THEN
776 0 : CPABORT("K points not handled yet for printing MO_CUBE")
777 : END IF
778 :
779 324 : dft_section => section_vals_get_subs_vals(input, "DFT")
780 324 : logger => cp_get_default_logger()
781 :
782 : ! Quick return if no print required
783 324 : IF (.NOT. BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
784 : "PRINT%MO_CUBES"), cp_p_file)) THEN
785 300 : CALL timestop(handle)
786 300 : RETURN
787 : END IF
788 :
789 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, &
790 : mos=mos, &
791 : blacs_env=blacs_env, &
792 : qs_kind_set=qs_kind_set, &
793 : pw_env=pw_env, &
794 : subsys=subsys, &
795 : para_env=para_env, &
796 : particle_set=particle_set, &
797 24 : dft_control=dft_control)
798 24 : CALL qs_subsys_get(subsys, particles=particles)
799 :
800 24 : nspins = dft_control%nspins
801 24 : rt_time_step = qs_env%sim_step
802 :
803 : ! Setup the grids needed to compute a wavefunction given a vector
804 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
805 24 : pw_pools=pw_pools)
806 24 : CALL auxbas_pw_pool%create_pw(wf_r)
807 24 : CALL auxbas_pw_pool%create_pw(wf_g)
808 24 : CALL auxbas_pw_pool%create_pw(density_r)
809 24 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
810 :
811 70 : DO ispin = 1, nspins
812 46 : CALL get_mo_set(mo_set=mos(ispin), homo=homo)
813 :
814 46 : nhomo = section_get_ival(dft_section, "PRINT%MO_CUBES%NHOMO")
815 46 : append_cube = section_get_lval(dft_section, "PRINT%MO_CUBES%APPEND")
816 46 : my_pos_cube = "REWIND"
817 46 : IF (append_cube) THEN
818 0 : my_pos_cube = "APPEND"
819 : END IF
820 46 : CALL section_vals_val_get(dft_section, "PRINT%MO_CUBES%HOMO_LIST", n_rep_val=n_rep)
821 46 : IF (n_rep > 0) THEN ! write the cubes of the list
822 0 : nlist = 0
823 0 : DO ir = 1, n_rep
824 0 : NULLIFY (list)
825 : CALL section_vals_val_get(dft_section, "PRINT%MO_CUBES%HOMO_LIST", i_rep_val=ir, &
826 0 : i_vals=list)
827 0 : IF (ASSOCIATED(list)) THEN
828 0 : CALL reallocate(list_index, 1, nlist + SIZE(list))
829 0 : DO i = 1, SIZE(list)
830 0 : list_index(i + nlist) = list(i)
831 : END DO
832 0 : nlist = nlist + SIZE(list)
833 : END IF
834 : END DO
835 : ELSE
836 :
837 46 : IF (nhomo == -1) nhomo = homo
838 46 : nlist = homo - MAX(1, homo - nhomo + 1) + 1
839 138 : ALLOCATE (list_index(nlist))
840 224 : DO i = 1, nlist
841 224 : list_index(i) = MAX(1, homo - nhomo + 1) + i - 1
842 : END DO
843 : END IF
844 224 : DO i = 1, nlist
845 178 : ivector = list_index(i)
846 : CALL get_qs_env(qs_env=qs_env, &
847 : atomic_kind_set=atomic_kind_set, &
848 : qs_kind_set=qs_kind_set, &
849 : cell=cell, &
850 : particle_set=particle_set, &
851 178 : pw_env=pw_env)
852 :
853 : ! density_r contains the density of the MOs
854 178 : CALL pw_zero(density_r)
855 178 : mo_coeff => mos_new(2*ispin - 1)!Real coeff
856 : CALL calculate_wavefunction(mo_coeff, ivector, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
857 178 : cell, dft_control, particle_set, pw_env)
858 : ! Adding the real part
859 178 : CALL pw_multiply(density_r, wf_r, wf_r, 1.0_dp)
860 :
861 178 : mo_coeff => mos_new(2*ispin) !Im coeff
862 : CALL calculate_wavefunction(mo_coeff, ivector, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
863 178 : cell, dft_control, particle_set, pw_env)
864 : ! Adding the im part
865 178 : CALL pw_multiply(density_r, wf_r, wf_r, 1.0_dp)
866 :
867 178 : WRITE (filename, '(a4,I5.5,a1,I1.1)') "WFN_", ivector, "_", ispin
868 178 : mpi_io = .TRUE.
869 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MO_CUBES", extension=".cube", &
870 : middle_name=TRIM(filename), file_position=my_pos_cube, log_filename=.FALSE., &
871 178 : mpi_io=mpi_io)
872 178 : WRITE (title, *) "DENSITY ", ivector, " spin ", ispin, " i.e. HOMO - ", ivector - homo
873 : CALL cp_pw_to_cube(density_r, unit_nr, title, particles=particles, &
874 178 : stride=section_get_ivals(dft_section, "PRINT%MO_CUBES%STRIDE"), mpi_io=mpi_io)
875 224 : CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MO_CUBES", mpi_io=mpi_io)
876 : END DO
877 162 : IF (ASSOCIATED(list_index)) DEALLOCATE (list_index)
878 : END DO
879 :
880 : ! Deallocate grids needed to compute wavefunctions
881 24 : CALL auxbas_pw_pool%give_back_pw(wf_r)
882 24 : CALL auxbas_pw_pool%give_back_pw(wf_g)
883 24 : CALL auxbas_pw_pool%give_back_pw(density_r)
884 :
885 24 : CALL timestop(handle)
886 :
887 324 : END SUBROUTINE write_rtp_mo_cubes
888 :
889 : ! **************************************************************************************************
890 : !> \brief Warn about unused sections of the print section - only implemented for some of the methods
891 : !> \param section Parent section
892 : !> \param subsection_name Name of the subsection which is not used even when explicitly included
893 : !> \param error_message Message to display as warning
894 : !> \author Stepan Marek (01.25)
895 : ! **************************************************************************************************
896 548 : SUBROUTINE warn_section_unused(section, subsection_name, error_message)
897 : TYPE(section_vals_type), POINTER :: section
898 : CHARACTER(len=*) :: subsection_name, error_message
899 :
900 : LOGICAL :: explicit
901 : TYPE(section_vals_type), POINTER :: target_section
902 :
903 274 : target_section => section_vals_get_subs_vals(section, subsection_name)
904 274 : CALL section_vals_get(target_section, explicit=explicit)
905 274 : IF (explicit) CPWARN(error_message)
906 274 : END SUBROUTINE warn_section_unused
907 : ! **************************************************************************************************
908 : !> \brief Attempt to read the moments from a previously written file
909 : !> \param moments_section Print key section for the moments trace
910 : !> \param orig_start Index from which the original calculation started (given by input file)
911 : !> \param current_start The current start index from which the calculation continues
912 : !> \param moments Moment trace, where the read entries are overwritten
913 : !> \param times Optional time trace, where again the read entries are overwritten
914 : !> \param mom_read Whether moments were actually read
915 : !> \author Stepan Marek (11.25)
916 : ! **************************************************************************************************
917 32 : SUBROUTINE read_moments(moments_section, orig_start, current_start, moments, times, mom_read)
918 : TYPE(section_vals_type), POINTER :: moments_section
919 : INTEGER :: orig_start, current_start
920 : COMPLEX(kind=dp), DIMENSION(:, :, :), POINTER :: moments
921 : REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: times
922 : LOGICAL, OPTIONAL :: mom_read
923 :
924 : CHARACTER(len=14), DIMENSION(4) :: file_extensions
925 : CHARACTER(len=default_path_length) :: save_name_im, save_name_re
926 : INTEGER :: i, k, max_n, n, n_spin, unit_im, unit_re
927 : REAL(kind=dp) :: time
928 : REAL(kind=dp), DIMENSION(3) :: moments_im, moments_re
929 : TYPE(cp_logger_type), POINTER :: logger
930 :
931 : ! Whether at least some moments where extracted from the files
932 :
933 32 : file_extensions(1) = "_SPIN_A_RE.dat"
934 32 : file_extensions(2) = "_SPIN_A_IM.dat"
935 32 : file_extensions(3) = "_SPIN_B_RE.dat"
936 32 : file_extensions(4) = "_SPIN_B_IM.dat"
937 :
938 32 : IF (PRESENT(mom_read)) mom_read = .FALSE.
939 :
940 32 : n_spin = SIZE(moments, 1)
941 :
942 32 : logger => cp_get_default_logger()
943 32 : max_n = current_start - orig_start + 1
944 :
945 76 : DO i = 1, n_spin
946 : ! Open the relevant files
947 : save_name_re = cp_print_key_generate_filename(logger, moments_section, &
948 44 : extension=file_extensions(2*i - 1), my_local=.FALSE.)
949 : save_name_im = cp_print_key_generate_filename(logger, moments_section, &
950 44 : extension=file_extensions(2*i), my_local=.FALSE.)
951 76 : IF (file_exists(save_name_re) .AND. file_exists(save_name_im)) THEN
952 : CALL open_file(save_name_re, file_status="OLD", file_form="FORMATTED", file_action="READ", &
953 0 : unit_number=unit_re)
954 : CALL open_file(save_name_im, file_status="OLD", file_form="FORMATTED", file_action="READ", &
955 0 : unit_number=unit_im)
956 : ! Skip the headers
957 0 : READ (unit_re, '()')
958 0 : READ (unit_im, '()')
959 0 : DO n = 1, max_n
960 : ! Read the data
961 0 : READ (unit_re, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3)') time, &
962 0 : moments_re(1), moments_re(2), moments_re(3)
963 0 : READ (unit_im, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3)') time, &
964 0 : moments_im(1), moments_im(2), moments_im(3)
965 0 : DO k = 1, 3
966 0 : moments(i, k, n) = CMPLX(moments_re(k), moments_im(k), kind=dp)
967 : END DO
968 0 : IF (PRESENT(times)) THEN
969 0 : times(n) = time/femtoseconds
970 : END IF
971 : END DO
972 0 : IF (PRESENT(mom_read)) mom_read = .TRUE.
973 0 : CALL close_file(unit_re)
974 0 : CALL close_file(unit_im)
975 : ELSE
976 948 : moments(i, :, 1:max_n) = CMPLX(0.0, 0.0, kind=dp)
977 44 : CPWARN("Could not read at least one moments file - missing trace is set to zero.")
978 : END IF
979 : END DO
980 32 : END SUBROUTINE read_moments
981 :
982 : ! **************************************************************************************************
983 : !> \brief Recalculates the field for index range given by orig_start and i_start,
984 : !> with times taken from the times array
985 : !> \param fields Array, where relevant indices are recalculated
986 : !> \param times Array of input times
987 : !> \param orig_start original start (in case of restarts)
988 : !> \param i_start current start
989 : !> \param dft_control DFT parameters
990 : !> \date 11.2025
991 : !> \author Stepan Marek
992 : ! **************************************************************************************************
993 18 : SUBROUTINE recalculate_fields(fields, times, orig_start, i_start, dft_control)
994 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: fields
995 : REAL(kind=dp), DIMENSION(:), POINTER :: times
996 : INTEGER :: orig_start, i_start
997 : TYPE(dft_control_type), POINTER :: dft_control
998 :
999 : INTEGER :: i, max_n
1000 : REAL(kind=dp), DIMENSION(3) :: field
1001 :
1002 18 : max_n = i_start - orig_start + 1
1003 :
1004 18 : IF (dft_control%rtp_control%apply_delta_pulse .OR. &
1005 : dft_control%rtp_control%apply_delta_pulse_mag) THEN
1006 90 : fields(:, 1:max_n) = 0.0_dp
1007 : ELSE
1008 0 : DO i = 1, max_n
1009 0 : CALL make_field(dft_control, field, i + orig_start - 1, times(i))
1010 0 : fields(:, i) = field(:)
1011 : END DO
1012 : END IF
1013 18 : END SUBROUTINE recalculate_fields
1014 :
1015 : END MODULE rt_propagation_utils
|
