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