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