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