(git:e7e05ae)
qs_rho_methods.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief methods of the rho structure (defined in qs_rho_types)
10 !> \par History
11 !> 08.2002 created [fawzi]
12 !> 08.2014 kpoints [JGH]
13 !> \author Fawzi Mohamed
14 ! **************************************************************************************************
15 MODULE qs_rho_methods
16  USE admm_types, ONLY: get_admm_env
17  USE atomic_kind_types, ONLY: atomic_kind_type
18  USE cp_control_types, ONLY: dft_control_type
19  USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
20  USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
21  dbcsr_deallocate_matrix_set
22  USE cp_log_handling, ONLY: cp_to_string
23  USE dbcsr_api, ONLY: &
24  dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_p_type, dbcsr_scale, dbcsr_set, dbcsr_type, &
25  dbcsr_type_antisymmetric, dbcsr_type_symmetric
26  USE kinds, ONLY: default_string_length,&
27  dp
28  USE kpoint_types, ONLY: get_kpoint_info,&
29  kpoint_type
30  USE lri_environment_methods, ONLY: calculate_lri_densities
31  USE lri_environment_types, ONLY: lri_density_type,&
32  lri_environment_type
33  USE message_passing, ONLY: mp_para_env_type
34  USE pw_env_types, ONLY: pw_env_get,&
35  pw_env_type
36  USE pw_methods, ONLY: pw_axpy,&
37  pw_copy,&
38  pw_scale,&
39  pw_zero
40  USE pw_pool_types, ONLY: pw_pool_type
41  USE pw_types, ONLY: pw_c1d_gs_type,&
42  pw_r3d_rs_type
43  USE qs_collocate_density, ONLY: calculate_drho_elec,&
44  calculate_rho_elec
45  USE qs_environment_types, ONLY: get_qs_env,&
46  qs_environment_type,&
47  set_qs_env
48  USE qs_kind_types, ONLY: qs_kind_type
49  USE qs_ks_types, ONLY: get_ks_env,&
50  qs_ks_env_type
51  USE qs_local_rho_types, ONLY: local_rho_type
52  USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
53  USE qs_oce_types, ONLY: oce_matrix_type
54  USE qs_rho_atom_methods, ONLY: calculate_rho_atom_coeff
55  USE qs_rho_atom_types, ONLY: rho_atom_type
56  USE qs_rho_types, ONLY: qs_rho_clear,&
57  qs_rho_get,&
58  qs_rho_set,&
59  qs_rho_type
60  USE ri_environment_methods, ONLY: calculate_ri_densities
61  USE task_list_types, ONLY: task_list_type
62 #include "./base/base_uses.f90"
63 
64  IMPLICIT NONE
65  PRIVATE
66 
67  LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
68  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho_methods'
69 
70  PUBLIC :: qs_rho_update_rho, qs_rho_update_tddfpt, &
71  qs_rho_rebuild, qs_rho_copy, qs_rho_scale_and_add
72  PUBLIC :: duplicate_rho_type, allocate_rho_ao_imag_from_real
73 
74 CONTAINS
75 
76 ! **************************************************************************************************
77 !> \brief rebuilds rho (if necessary allocating and initializing it)
78 !> \param rho the rho type to rebuild (defaults to qs_env%rho)
79 !> \param qs_env the environment to which rho belongs
80 !> \param rebuild_ao if it is necessary to rebuild rho_ao. Defaults to true.
81 !> \param rebuild_grids if it in necessary to rebuild rho_r and rho_g.
82 !> Defaults to false.
83 !> \param admm (use aux_fit basis)
84 !> \param pw_env_external external plane wave environment
85 !> \par History
86 !> 11.2002 created replacing qs_rho_create and qs_env_rebuild_rho[fawzi]
87 !> \author Fawzi Mohamed
88 !> \note
89 !> needs updated pw pools, s, s_mstruct and h in qs_env.
90 !> The use of p to keep the structure of h (needed for the forces)
91 !> is ugly and should be removed.
92 !> Change so that it does not allocate a subcomponent if it is not
93 !> associated and not requested?
94 ! **************************************************************************************************
95  SUBROUTINE qs_rho_rebuild(rho, qs_env, rebuild_ao, rebuild_grids, admm, pw_env_external)
96  TYPE(qs_rho_type), INTENT(INOUT) :: rho
97  TYPE(qs_environment_type), POINTER :: qs_env
98  LOGICAL, INTENT(in), OPTIONAL :: rebuild_ao, rebuild_grids, admm
99  TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
100 
101  CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_rho_rebuild'
102 
103  CHARACTER(LEN=default_string_length) :: headline
104  INTEGER :: handle, i, ic, j, nimg, nspins
105  LOGICAL :: do_kpoints, my_admm, my_rebuild_ao, &
106  my_rebuild_grids, rho_ao_is_complex
107  REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r
108  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, rho_ao_im_kp, rho_ao_kp
109  TYPE(dbcsr_type), POINTER :: refmatrix, tmatrix
110  TYPE(dft_control_type), POINTER :: dft_control
111  TYPE(kpoint_type), POINTER :: kpoints
112  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
113  POINTER :: sab_orb
114  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, tau_g
115  TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER :: drho_g
116  TYPE(pw_env_type), POINTER :: pw_env
117  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
118  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau_r
119  TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER :: drho_r
120  TYPE(pw_r3d_rs_type), POINTER :: rho_r_sccs
121 
122  CALL timeset(routinen, handle)
123 
124  NULLIFY (pw_env, auxbas_pw_pool, matrix_s_kp, dft_control)
125  NULLIFY (tot_rho_r, rho_ao_kp, rho_r, rho_g, drho_r, drho_g, tau_r, tau_g, rho_ao_im_kp)
126  NULLIFY (rho_r_sccs)
127  NULLIFY (sab_orb)
128  my_rebuild_ao = .true.
129  my_rebuild_grids = .true.
130  my_admm = .false.
131  IF (PRESENT(rebuild_ao)) my_rebuild_ao = rebuild_ao
132  IF (PRESENT(rebuild_grids)) my_rebuild_grids = rebuild_grids
133  IF (PRESENT(admm)) my_admm = admm
134 
135  CALL get_qs_env(qs_env, &
136  kpoints=kpoints, &
137  do_kpoints=do_kpoints, &
138  pw_env=pw_env, &
139  dft_control=dft_control)
140  IF (PRESENT(pw_env_external)) &
141  pw_env => pw_env_external
142 
143  nimg = dft_control%nimages
144 
145  IF (my_admm) THEN
146  CALL get_admm_env(qs_env%admm_env, sab_aux_fit=sab_orb, matrix_s_aux_fit_kp=matrix_s_kp)
147  ELSE
148  CALL get_qs_env(qs_env, matrix_s_kp=matrix_s_kp)
149 
150  IF (do_kpoints) THEN
151  CALL get_kpoint_info(kpoints, sab_nl=sab_orb)
152  ELSE
153  CALL get_qs_env(qs_env, sab_orb=sab_orb)
154  END IF
155  END IF
156  refmatrix => matrix_s_kp(1, 1)%matrix
157 
158  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
159  nspins = dft_control%nspins
160 
161  CALL qs_rho_get(rho, &
162  tot_rho_r=tot_rho_r, &
163  rho_ao_kp=rho_ao_kp, &
164  rho_ao_im_kp=rho_ao_im_kp, &
165  rho_r=rho_r, &
166  rho_g=rho_g, &
167  drho_r=drho_r, &
168  drho_g=drho_g, &
169  tau_r=tau_r, &
170  tau_g=tau_g, &
171  rho_r_sccs=rho_r_sccs, &
172  complex_rho_ao=rho_ao_is_complex)
173 
174  IF (.NOT. ASSOCIATED(tot_rho_r)) THEN
175  ALLOCATE (tot_rho_r(nspins))
176  tot_rho_r = 0.0_dp
177  CALL qs_rho_set(rho, tot_rho_r=tot_rho_r)
178  END IF
179 
180  ! rho_ao
181  IF (my_rebuild_ao .OR. (.NOT. ASSOCIATED(rho_ao_kp))) THEN
182  IF (ASSOCIATED(rho_ao_kp)) &
183  CALL dbcsr_deallocate_matrix_set(rho_ao_kp)
184  ! Create a new density matrix set
185  CALL dbcsr_allocate_matrix_set(rho_ao_kp, nspins, nimg)
186  CALL qs_rho_set(rho, rho_ao_kp=rho_ao_kp)
187  DO i = 1, nspins
188  DO ic = 1, nimg
189  IF (nspins > 1) THEN
190  IF (i == 1) THEN
191  headline = "DENSITY MATRIX FOR ALPHA SPIN"
192  ELSE
193  headline = "DENSITY MATRIX FOR BETA SPIN"
194  END IF
195  ELSE
196  headline = "DENSITY MATRIX"
197  END IF
198  ALLOCATE (rho_ao_kp(i, ic)%matrix)
199  tmatrix => rho_ao_kp(i, ic)%matrix
200  CALL dbcsr_create(matrix=tmatrix, template=refmatrix, name=trim(headline), &
201  matrix_type=dbcsr_type_symmetric, nze=0)
202  CALL cp_dbcsr_alloc_block_from_nbl(tmatrix, sab_orb)
203  CALL dbcsr_set(tmatrix, 0.0_dp)
204  END DO
205  END DO
206  IF (rho_ao_is_complex) THEN
207  IF (ASSOCIATED(rho_ao_im_kp)) THEN
208  CALL dbcsr_deallocate_matrix_set(rho_ao_im_kp)
209  END IF
210  CALL dbcsr_allocate_matrix_set(rho_ao_im_kp, nspins, nimg)
211  CALL qs_rho_set(rho, rho_ao_im_kp=rho_ao_im_kp)
212  DO i = 1, nspins
213  DO ic = 1, nimg
214  IF (nspins > 1) THEN
215  IF (i == 1) THEN
216  headline = "IMAGINARY PART OF DENSITY MATRIX FOR ALPHA SPIN"
217  ELSE
218  headline = "IMAGINARY PART OF DENSITY MATRIX FOR BETA SPIN"
219  END IF
220  ELSE
221  headline = "IMAGINARY PART OF DENSITY MATRIX"
222  END IF
223  ALLOCATE (rho_ao_im_kp(i, ic)%matrix)
224  tmatrix => rho_ao_im_kp(i, ic)%matrix
225  CALL dbcsr_create(matrix=tmatrix, template=refmatrix, name=trim(headline), &
226  matrix_type=dbcsr_type_antisymmetric, nze=0)
227  CALL cp_dbcsr_alloc_block_from_nbl(tmatrix, sab_orb)
228  CALL dbcsr_set(tmatrix, 0.0_dp)
229  END DO
230  END DO
231  END IF
232  END IF
233 
234  ! rho_r
235  IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(rho_r)) THEN
236  IF (ASSOCIATED(rho_r)) THEN
237  DO i = 1, SIZE(rho_r)
238  CALL rho_r(i)%release()
239  END DO
240  DEALLOCATE (rho_r)
241  END IF
242  ALLOCATE (rho_r(nspins))
243  CALL qs_rho_set(rho, rho_r=rho_r)
244  DO i = 1, nspins
245  CALL auxbas_pw_pool%create_pw(rho_r(i))
246  END DO
247  END IF
248 
249  ! rho_g
250  IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(rho_g)) THEN
251  IF (ASSOCIATED(rho_g)) THEN
252  DO i = 1, SIZE(rho_g)
253  CALL rho_g(i)%release()
254  END DO
255  DEALLOCATE (rho_g)
256  END IF
257  ALLOCATE (rho_g(nspins))
258  CALL qs_rho_set(rho, rho_g=rho_g)
259  DO i = 1, nspins
260  CALL auxbas_pw_pool%create_pw(rho_g(i))
261  END DO
262  END IF
263 
264  ! SCCS
265  IF (dft_control%do_sccs) THEN
266  IF (my_rebuild_grids .OR. (.NOT. ASSOCIATED(rho_r_sccs))) THEN
267  IF (ASSOCIATED(rho_r_sccs)) THEN
268  CALL rho_r_sccs%release()
269  DEALLOCATE (rho_r_sccs)
270  END IF
271  ALLOCATE (rho_r_sccs)
272  CALL qs_rho_set(rho, rho_r_sccs=rho_r_sccs)
273  CALL auxbas_pw_pool%create_pw(rho_r_sccs)
274  CALL pw_zero(rho_r_sccs)
275  END IF
276  END IF
277 
278  ! allocate drho_r and drho_g if xc_deriv_collocate
279  IF (dft_control%drho_by_collocation) THEN
280  ! drho_r
281  IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(drho_r)) THEN
282  IF (ASSOCIATED(drho_r)) THEN
283  DO j = 1, SIZE(drho_r, 2)
284  DO i = 1, SIZE(drho_r, 1)
285  CALL drho_r(i, j)%release()
286  END DO
287  END DO
288  DEALLOCATE (drho_r)
289  END IF
290  ALLOCATE (drho_r(3, nspins))
291  CALL qs_rho_set(rho, drho_r=drho_r)
292  DO j = 1, nspins
293  DO i = 1, 3
294  CALL auxbas_pw_pool%create_pw(drho_r(i, j))
295  END DO
296  END DO
297  END IF
298  ! drho_g
299  IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(drho_g)) THEN
300  IF (ASSOCIATED(drho_g)) THEN
301  DO j = 1, SIZE(drho_g, 2)
302  DO i = 1, SIZE(drho_r, 1)
303  CALL drho_g(i, j)%release()
304  END DO
305  END DO
306  DEALLOCATE (drho_g)
307  END IF
308  ALLOCATE (drho_g(3, nspins))
309  CALL qs_rho_set(rho, drho_g=drho_g)
310  DO j = 1, nspins
311  DO i = 1, 3
312  CALL auxbas_pw_pool%create_pw(drho_g(i, j))
313  END DO
314  END DO
315  END IF
316  END IF
317 
318  ! allocate tau_r and tau_g if use_kinetic_energy_density
319  IF (dft_control%use_kinetic_energy_density) THEN
320  ! tau_r
321  IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(tau_r)) THEN
322  IF (ASSOCIATED(tau_r)) THEN
323  DO i = 1, SIZE(tau_r)
324  CALL tau_r(i)%release()
325  END DO
326  DEALLOCATE (tau_r)
327  END IF
328  ALLOCATE (tau_r(nspins))
329  CALL qs_rho_set(rho, tau_r=tau_r)
330  DO i = 1, nspins
331  CALL auxbas_pw_pool%create_pw(tau_r(i))
332  END DO
333  END IF
334 
335  ! tau_g
336  IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(tau_g)) THEN
337  IF (ASSOCIATED(tau_g)) THEN
338  DO i = 1, SIZE(tau_g)
339  CALL tau_g(i)%release()
340  END DO
341  DEALLOCATE (tau_g)
342  END IF
343  ALLOCATE (tau_g(nspins))
344  CALL qs_rho_set(rho, tau_g=tau_g)
345  DO i = 1, nspins
346  CALL auxbas_pw_pool%create_pw(tau_g(i))
347  END DO
348  END IF
349  END IF ! use_kinetic_energy_density
350 
351  CALL timestop(handle)
352 
353  END SUBROUTINE qs_rho_rebuild
354 
355 ! **************************************************************************************************
356 !> \brief updates rho_r and rho_g to the rho%rho_ao.
357 !> if use_kinetic_energy_density also computes tau_r and tau_g
358 !> this works for all ground state and ground state response methods
359 !> \param rho_struct the rho structure that should be updated
360 !> \param qs_env the qs_env rho_struct refers to
361 !> the integrated charge in r space
362 !> \param rho_xc_external ...
363 !> \param local_rho_set ...
364 !> \param task_list_external external task list
365 !> \param task_list_external_soft external task list (soft_version)
366 !> \param pw_env_external external plane wave environment
367 !> \param para_env_external external MPI environment
368 !> \par History
369 !> 08.2002 created [fawzi]
370 !> \author Fawzi Mohamed
371 ! **************************************************************************************************
372  SUBROUTINE qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, &
373  task_list_external, task_list_external_soft, &
374  pw_env_external, para_env_external)
375  TYPE(qs_rho_type), INTENT(INOUT) :: rho_struct
376  TYPE(qs_environment_type), POINTER :: qs_env
377  TYPE(qs_rho_type), OPTIONAL, POINTER :: rho_xc_external
378  TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set
379  TYPE(task_list_type), OPTIONAL, POINTER :: task_list_external, &
380  task_list_external_soft
381  TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
382  TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_external
383 
384  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
385  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
386  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
387  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
388  TYPE(dft_control_type), POINTER :: dft_control
389  TYPE(kpoint_type), POINTER :: kpoints
390  TYPE(lri_density_type), POINTER :: lri_density
391  TYPE(lri_environment_type), POINTER :: lri_env
392  TYPE(mp_para_env_type), POINTER :: para_env
393  TYPE(qs_ks_env_type), POINTER :: ks_env
394 
395  CALL get_qs_env(qs_env, dft_control=dft_control, &
396  atomic_kind_set=atomic_kind_set, &
397  para_env=para_env)
398  IF (PRESENT(para_env_external)) para_env => para_env_external
399 
400  IF (dft_control%qs_control%semi_empirical .OR. &
401  dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
402 
403  CALL qs_rho_set(rho_struct, rho_r_valid=.false., rho_g_valid=.false.)
404 
405  ELSEIF (dft_control%qs_control%lrigpw) THEN
406  cpassert(.NOT. dft_control%use_kinetic_energy_density)
407  cpassert(.NOT. dft_control%drho_by_collocation)
408  CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
409  CALL get_qs_env(qs_env, ks_env=ks_env)
410  CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
411  CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
412  CALL get_qs_env(qs_env, lri_env=lri_env, lri_density=lri_density)
413  CALL calculate_lri_densities(lri_env, lri_density, qs_env, rho_ao_kp, cell_to_index, &
414  lri_rho_struct=rho_struct, &
415  atomic_kind_set=atomic_kind_set, &
416  para_env=para_env, &
417  response_density=.false.)
418  CALL set_qs_env(qs_env, lri_density=lri_density)
419  CALL qs_rho_set(rho_struct, rho_r_valid=.true., rho_g_valid=.true.)
420 
421  ELSEIF (dft_control%qs_control%rigpw) THEN
422  cpassert(.NOT. dft_control%use_kinetic_energy_density)
423  cpassert(.NOT. dft_control%drho_by_collocation)
424  CALL get_qs_env(qs_env, lri_env=lri_env)
425  CALL qs_rho_get(rho_struct, rho_ao=rho_ao)
426  CALL calculate_ri_densities(lri_env, qs_env, rho_ao, &
427  lri_rho_struct=rho_struct, &
428  atomic_kind_set=atomic_kind_set, &
429  para_env=para_env)
430  CALL qs_rho_set(rho_struct, rho_r_valid=.true., rho_g_valid=.true.)
431 
432  ELSE
433  CALL qs_rho_update_rho_low(rho_struct=rho_struct, qs_env=qs_env, &
434  rho_xc_external=rho_xc_external, &
435  local_rho_set=local_rho_set, &
436  task_list_external=task_list_external, &
437  task_list_external_soft=task_list_external_soft, &
438  pw_env_external=pw_env_external, &
439  para_env_external=para_env_external)
440 
441  END IF
442 
443  END SUBROUTINE qs_rho_update_rho
444 
445 ! **************************************************************************************************
446 !> \brief updates rho_r and rho_g to the rho%rho_ao.
447 !> if use_kinetic_energy_density also computes tau_r and tau_g
448 !> \param rho_struct the rho structure that should be updated
449 !> \param qs_env the qs_env rho_struct refers to
450 !> the integrated charge in r space
451 !> \param rho_xc_external rho structure for GAPW_XC
452 !> \param local_rho_set ...
453 !> \param pw_env_external external plane wave environment
454 !> \param task_list_external external task list (use for default and GAPW)
455 !> \param task_list_external_soft external task list (soft density for GAPW_XC)
456 !> \param para_env_external ...
457 !> \par History
458 !> 08.2002 created [fawzi]
459 !> \author Fawzi Mohamed
460 ! **************************************************************************************************
461  SUBROUTINE qs_rho_update_rho_low(rho_struct, qs_env, rho_xc_external, &
462  local_rho_set, pw_env_external, &
463  task_list_external, task_list_external_soft, &
464  para_env_external)
465  TYPE(qs_rho_type), INTENT(INOUT) :: rho_struct
466  TYPE(qs_environment_type), POINTER :: qs_env
467  TYPE(qs_rho_type), OPTIONAL, POINTER :: rho_xc_external
468  TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set
469  TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
470  TYPE(task_list_type), OPTIONAL, POINTER :: task_list_external, &
471  task_list_external_soft
472  TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_external
473 
474  CHARACTER(len=*), PARAMETER :: routinen = 'qs_rho_update_rho_low'
475 
476  INTEGER :: handle, img, ispin, nimg, nspins
477  LOGICAL :: gapw, gapw_xc
478  REAL(kind=dp) :: dum
479  REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r, tot_rho_r_xc
480  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
481  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
482  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp, rho_xc_ao
483  TYPE(dft_control_type), POINTER :: dft_control
484  TYPE(mp_para_env_type), POINTER :: para_env
485  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
486  POINTER :: sab
487  TYPE(oce_matrix_type), POINTER :: oce
488  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rho_xc_g, tau_g, tau_xc_g
489  TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER :: drho_g, drho_xc_g
490  TYPE(pw_env_type), POINTER :: pw_env
491  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, rho_xc_r, tau_r, tau_xc_r
492  TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER :: drho_r, drho_xc_r
493  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
494  TYPE(qs_ks_env_type), POINTER :: ks_env
495  TYPE(qs_rho_type), POINTER :: rho_xc
496  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
497  TYPE(task_list_type), POINTER :: task_list
498 
499  CALL timeset(routinen, handle)
500 
501  NULLIFY (dft_control, rho_xc, ks_env, rho_ao, rho_r, rho_g, drho_r, drho_g, tau_r, tau_g)
502  NULLIFY (rho_xc_ao, rho_xc_g, rho_xc_r, drho_xc_g, tau_xc_r, tau_xc_g, tot_rho_r, tot_rho_r_xc)
503  NULLIFY (para_env, pw_env, atomic_kind_set)
504 
505  CALL get_qs_env(qs_env, &
506  ks_env=ks_env, &
507  dft_control=dft_control, &
508  atomic_kind_set=atomic_kind_set)
509 
510  CALL qs_rho_get(rho_struct, &
511  rho_r=rho_r, &
512  rho_g=rho_g, &
513  tot_rho_r=tot_rho_r, &
514  drho_r=drho_r, &
515  drho_g=drho_g, &
516  tau_r=tau_r, &
517  tau_g=tau_g)
518 
519  CALL get_qs_env(qs_env, task_list=task_list, &
520  para_env=para_env, pw_env=pw_env)
521  IF (PRESENT(pw_env_external)) pw_env => pw_env_external
522  IF (PRESENT(task_list_external)) task_list => task_list_external
523  IF (PRESENT(para_env_external)) para_env => para_env_external
524 
525  nspins = dft_control%nspins
526  nimg = dft_control%nimages
527  gapw = dft_control%qs_control%gapw
528  gapw_xc = dft_control%qs_control%gapw_xc
529 
530  CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
531  DO ispin = 1, nspins
532  rho_ao => rho_ao_kp(ispin, :)
533  CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
534  rho=rho_r(ispin), &
535  rho_gspace=rho_g(ispin), &
536  total_rho=tot_rho_r(ispin), &
537  ks_env=ks_env, soft_valid=gapw, &
538  task_list_external=task_list_external, &
539  pw_env_external=pw_env_external)
540  END DO
541  CALL qs_rho_set(rho_struct, rho_r_valid=.true., rho_g_valid=.true.)
542 
543  IF (gapw_xc) THEN
544  IF (PRESENT(rho_xc_external)) THEN
545  rho_xc => rho_xc_external
546  ELSE
547  CALL get_qs_env(qs_env=qs_env, rho_xc=rho_xc)
548  END IF
549  CALL qs_rho_get(rho_xc, &
550  rho_ao_kp=rho_xc_ao, &
551  rho_r=rho_xc_r, &
552  rho_g=rho_xc_g, &
553  tot_rho_r=tot_rho_r_xc)
554  ! copy rho_ao into rho_xc_ao
555  DO ispin = 1, nspins
556  DO img = 1, nimg
557  CALL dbcsr_copy(rho_xc_ao(ispin, img)%matrix, rho_ao_kp(ispin, img)%matrix)
558  END DO
559  END DO
560  DO ispin = 1, nspins
561  rho_ao => rho_xc_ao(ispin, :)
562  CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
563  rho=rho_xc_r(ispin), &
564  rho_gspace=rho_xc_g(ispin), &
565  total_rho=tot_rho_r_xc(ispin), &
566  ks_env=ks_env, soft_valid=gapw_xc, &
567  task_list_external=task_list_external_soft, &
568  pw_env_external=pw_env_external)
569  END DO
570  CALL qs_rho_set(rho_xc, rho_r_valid=.true., rho_g_valid=.true.)
571  END IF
572 
573  ! GAPW o GAPW_XC require the calculation of hard and soft local densities
574  IF (gapw .OR. gapw_xc) THEN
575  CALL get_qs_env(qs_env=qs_env, &
576  rho_atom_set=rho_atom_set, &
577  qs_kind_set=qs_kind_set, &
578  oce=oce, sab_orb=sab)
579  IF (PRESENT(local_rho_set)) rho_atom_set => local_rho_set%rho_atom_set
580  cpassert(ASSOCIATED(rho_atom_set))
581  CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
582  CALL calculate_rho_atom_coeff(qs_env, rho_ao_kp, rho_atom_set, qs_kind_set, oce, sab, para_env)
583  END IF
584 
585  IF (.NOT. gapw_xc) THEN
586  ! if needed compute also the gradient of the density
587  IF (dft_control%drho_by_collocation) THEN
588  CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
589  cpassert(.NOT. PRESENT(task_list_external))
590  DO ispin = 1, nspins
591  rho_ao => rho_ao_kp(ispin, :)
592  CALL calculate_drho_elec(matrix_p_kp=rho_ao, &
593  drho=drho_r(:, ispin), &
594  drho_gspace=drho_g(:, ispin), &
595  qs_env=qs_env, soft_valid=gapw)
596  END DO
597  CALL qs_rho_set(rho_struct, drho_r_valid=.true., drho_g_valid=.true.)
598  END IF
599  ! if needed compute also the kinetic energy density
600  IF (dft_control%use_kinetic_energy_density) THEN
601  CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
602  DO ispin = 1, nspins
603  rho_ao => rho_ao_kp(ispin, :)
604  CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
605  rho=tau_r(ispin), &
606  rho_gspace=tau_g(ispin), &
607  total_rho=dum, & ! presumably not meaningful
608  ks_env=ks_env, soft_valid=gapw, &
609  compute_tau=.true., &
610  task_list_external=task_list_external, &
611  pw_env_external=pw_env_external)
612  END DO
613  CALL qs_rho_set(rho_struct, tau_r_valid=.true., tau_g_valid=.true.)
614  END IF
615  ELSE
616  CALL qs_rho_get(rho_xc, &
617  drho_r=drho_xc_r, &
618  drho_g=drho_xc_g, &
619  tau_r=tau_xc_r, &
620  tau_g=tau_xc_g)
621  ! if needed compute also the gradient of the density
622  IF (dft_control%drho_by_collocation) THEN
623  cpassert(.NOT. PRESENT(task_list_external))
624  DO ispin = 1, nspins
625  rho_ao => rho_xc_ao(ispin, :)
626  CALL calculate_drho_elec(matrix_p_kp=rho_ao, &
627  drho=drho_xc_r(:, ispin), &
628  drho_gspace=drho_xc_g(:, ispin), &
629  qs_env=qs_env, soft_valid=gapw_xc)
630  END DO
631  CALL qs_rho_set(rho_xc, drho_r_valid=.true., drho_g_valid=.true.)
632  END IF
633  ! if needed compute also the kinetic energy density
634  IF (dft_control%use_kinetic_energy_density) THEN
635  DO ispin = 1, nspins
636  rho_ao => rho_xc_ao(ispin, :)
637  CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
638  rho=tau_xc_r(ispin), &
639  rho_gspace=tau_xc_g(ispin), &
640  ks_env=ks_env, soft_valid=gapw_xc, &
641  compute_tau=.true., &
642  task_list_external=task_list_external_soft, &
643  pw_env_external=pw_env_external)
644  END DO
645  CALL qs_rho_set(rho_xc, tau_r_valid=.true., tau_g_valid=.true.)
646  END IF
647  END IF
648 
649  CALL timestop(handle)
650 
651  END SUBROUTINE qs_rho_update_rho_low
652 
653 ! **************************************************************************************************
654 !> \brief updates rho_r and rho_g to the rho%rho_ao.
655 !> if use_kinetic_energy_density also computes tau_r and tau_g
656 !> \param rho_struct the rho structure that should be updated
657 !> \param qs_env the qs_env rho_struct refers to
658 !> the integrated charge in r space
659 !> \param pw_env_external external plane wave environment
660 !> \param task_list_external external task list
661 !> \param para_env_external ...
662 !> \param tddfpt_lri_env ...
663 !> \param tddfpt_lri_density ...
664 !> \par History
665 !> 08.2002 created [fawzi]
666 !> \author Fawzi Mohamed
667 ! **************************************************************************************************
668  SUBROUTINE qs_rho_update_tddfpt(rho_struct, qs_env, pw_env_external, task_list_external, &
669  para_env_external, tddfpt_lri_env, tddfpt_lri_density)
670  TYPE(qs_rho_type), INTENT(INOUT) :: rho_struct
671  TYPE(qs_environment_type), POINTER :: qs_env
672  TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
673  TYPE(task_list_type), OPTIONAL, POINTER :: task_list_external
674  TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_external
675  TYPE(lri_environment_type), OPTIONAL, POINTER :: tddfpt_lri_env
676  TYPE(lri_density_type), OPTIONAL, POINTER :: tddfpt_lri_density
677 
678  CHARACTER(len=*), PARAMETER :: routinen = 'qs_rho_update_tddfpt'
679 
680  INTEGER :: handle, ispin, nspins
681  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
682  LOGICAL :: lri_response
683  REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r
684  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
685  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
686  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
687  TYPE(dft_control_type), POINTER :: dft_control
688  TYPE(kpoint_type), POINTER :: kpoints
689  TYPE(mp_para_env_type), POINTER :: para_env
690  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
691  TYPE(pw_env_type), POINTER :: pw_env
692  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
693  TYPE(qs_ks_env_type), POINTER :: ks_env
694  TYPE(task_list_type), POINTER :: task_list
695 
696  CALL timeset(routinen, handle)
697 
698  CALL get_qs_env(qs_env, &
699  ks_env=ks_env, &
700  dft_control=dft_control, &
701  atomic_kind_set=atomic_kind_set, &
702  task_list=task_list, &
703  para_env=para_env, &
704  pw_env=pw_env)
705  IF (PRESENT(pw_env_external)) pw_env => pw_env_external
706  IF (PRESENT(task_list_external)) task_list => task_list_external
707  IF (PRESENT(para_env_external)) para_env => para_env_external
708 
709  CALL qs_rho_get(rho_struct, &
710  rho_r=rho_r, &
711  rho_g=rho_g, &
712  tot_rho_r=tot_rho_r)
713 
714  nspins = dft_control%nspins
715 
716  lri_response = PRESENT(tddfpt_lri_env)
717  IF (lri_response) THEN
718  cpassert(PRESENT(tddfpt_lri_density))
719  END IF
720 
721  cpassert(.NOT. dft_control%drho_by_collocation)
722  cpassert(.NOT. dft_control%use_kinetic_energy_density)
723  cpassert(.NOT. dft_control%qs_control%gapw)
724  cpassert(.NOT. dft_control%qs_control%gapw_xc)
725 
726  IF (lri_response) THEN
727  CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
728  CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
729  CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
730  CALL calculate_lri_densities(tddfpt_lri_env, tddfpt_lri_density, qs_env, rho_ao_kp, cell_to_index, &
731  lri_rho_struct=rho_struct, &
732  atomic_kind_set=atomic_kind_set, &
733  para_env=para_env, &
734  response_density=lri_response)
735  CALL qs_rho_set(rho_struct, rho_r_valid=.true., rho_g_valid=.true.)
736  ELSE
737  CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
738  DO ispin = 1, nspins
739  rho_ao => rho_ao_kp(ispin, :)
740  CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
741  rho=rho_r(ispin), &
742  rho_gspace=rho_g(ispin), &
743  total_rho=tot_rho_r(ispin), &
744  ks_env=ks_env, &
745  task_list_external=task_list_external, &
746  pw_env_external=pw_env_external)
747  END DO
748  CALL qs_rho_set(rho_struct, rho_r_valid=.true., rho_g_valid=.true.)
749  END IF
750 
751  CALL timestop(handle)
752 
753  END SUBROUTINE qs_rho_update_tddfpt
754 
755 ! **************************************************************************************************
756 !> \brief Allocate a density structure and fill it with data from an input structure
757 !> SIZE(rho_input) == mspin == 1 direct copy
758 !> SIZE(rho_input) == mspin == 2 direct copy of alpha and beta spin
759 !> SIZE(rho_input) == 1 AND mspin == 2 copy rho/2 into alpha and beta spin
760 !> \param rho_input ...
761 !> \param rho_output ...
762 !> \param auxbas_pw_pool ...
763 !> \param mspin ...
764 ! **************************************************************************************************
765  SUBROUTINE qs_rho_copy(rho_input, rho_output, auxbas_pw_pool, mspin)
766 
767  TYPE(qs_rho_type), INTENT(IN) :: rho_input
768  TYPE(qs_rho_type), INTENT(INOUT) :: rho_output
769  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
770  INTEGER, INTENT(IN) :: mspin
771 
772  CHARACTER(len=*), PARAMETER :: routinen = 'qs_rho_copy'
773 
774  INTEGER :: handle, i, j, nspins
775  LOGICAL :: complex_rho_ao, drho_g_valid_in, drho_r_valid_in, rho_g_valid_in, rho_r_valid_in, &
776  soft_valid_in, tau_g_valid_in, tau_r_valid_in
777  REAL(kind=dp) :: ospin
778  REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_g_in, tot_rho_g_out, &
779  tot_rho_r_in, tot_rho_r_out
780  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_im_in, rho_ao_im_out, rho_ao_in, &
781  rho_ao_out
782  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp_in
783  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_in, rho_g_out, tau_g_in, tau_g_out
784  TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER :: drho_g_in, drho_g_out
785  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_in, rho_r_out, tau_r_in, tau_r_out
786  TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER :: drho_r_in, drho_r_out
787  TYPE(pw_r3d_rs_type), POINTER :: rho_r_sccs_in, rho_r_sccs_out
788 
789  CALL timeset(routinen, handle)
790 
791  cpassert(mspin == 1 .OR. mspin == 2)
792  ospin = 1._dp/real(mspin, kind=dp)
793 
794  CALL qs_rho_clear(rho_output)
795 
796  NULLIFY (rho_ao_in, rho_ao_kp_in, rho_ao_im_in, rho_r_in, rho_g_in, drho_r_in, &
797  drho_g_in, tau_r_in, tau_g_in, tot_rho_r_in, tot_rho_g_in, rho_r_sccs_in)
798 
799  CALL qs_rho_get(rho_input, &
800  rho_ao=rho_ao_in, &
801  rho_ao_kp=rho_ao_kp_in, &
802  rho_ao_im=rho_ao_im_in, &
803  rho_r=rho_r_in, &
804  rho_g=rho_g_in, &
805  drho_r=drho_r_in, &
806  drho_g=drho_g_in, &
807  tau_r=tau_r_in, &
808  tau_g=tau_g_in, &
809  tot_rho_r=tot_rho_r_in, &
810  tot_rho_g=tot_rho_g_in, &
811  rho_g_valid=rho_g_valid_in, &
812  rho_r_valid=rho_r_valid_in, &
813  drho_g_valid=drho_g_valid_in, &
814  drho_r_valid=drho_r_valid_in, &
815  tau_r_valid=tau_r_valid_in, &
816  tau_g_valid=tau_g_valid_in, &
817  rho_r_sccs=rho_r_sccs_in, &
818  soft_valid=soft_valid_in, &
819  complex_rho_ao=complex_rho_ao)
820 
821  NULLIFY (rho_ao_out, rho_ao_im_out, rho_r_out, rho_g_out, drho_r_out, &
822  drho_g_out, tau_r_out, tau_g_out, tot_rho_r_out, tot_rho_g_out, rho_r_sccs_out)
823  ! rho_ao
824  IF (ASSOCIATED(rho_ao_in)) THEN
825  nspins = SIZE(rho_ao_in)
826  cpassert(mspin >= nspins)
827  CALL dbcsr_allocate_matrix_set(rho_ao_out, mspin)
828  CALL qs_rho_set(rho_output, rho_ao=rho_ao_out)
829  IF (mspin > nspins) THEN
830  DO i = 1, mspin
831  ALLOCATE (rho_ao_out(i)%matrix)
832  CALL dbcsr_copy(rho_ao_out(i)%matrix, rho_ao_in(1)%matrix, name="RHO copy")
833  CALL dbcsr_scale(rho_ao_out(i)%matrix, ospin)
834  END DO
835  ELSE
836  DO i = 1, nspins
837  ALLOCATE (rho_ao_out(i)%matrix)
838  CALL dbcsr_copy(rho_ao_out(i)%matrix, rho_ao_in(i)%matrix, name="RHO copy")
839  END DO
840  END IF
841  END IF
842 
843  ! rho_ao_kp
844  ! only for non-kp, we could probably just copy this pointer, should work also for non-kp?
845  !IF (ASSOCIATED(rho_ao_kp_in)) THEN
846  ! CPABORT("Copy not available")
847  !END IF
848 
849  ! rho_ao_im
850  IF (ASSOCIATED(rho_ao_im_in)) THEN
851  nspins = SIZE(rho_ao_im_in)
852  cpassert(mspin >= nspins)
853  CALL dbcsr_allocate_matrix_set(rho_ao_im_out, mspin)
854  CALL qs_rho_set(rho_output, rho_ao_im=rho_ao_im_out)
855  IF (mspin > nspins) THEN
856  DO i = 1, mspin
857  ALLOCATE (rho_ao_im_out(i)%matrix)
858  CALL dbcsr_copy(rho_ao_im_out(i)%matrix, rho_ao_im_in(1)%matrix, name="RHO copy")
859  CALL dbcsr_scale(rho_ao_im_out(i)%matrix, ospin)
860  END DO
861  ELSE
862  DO i = 1, nspins
863  ALLOCATE (rho_ao_im_out(i)%matrix)
864  CALL dbcsr_copy(rho_ao_im_out(i)%matrix, rho_ao_im_in(i)%matrix, name="RHO copy")
865  END DO
866  END IF
867  END IF
868 
869  ! rho_r
870  IF (ASSOCIATED(rho_r_in)) THEN
871  nspins = SIZE(rho_r_in)
872  cpassert(mspin >= nspins)
873  ALLOCATE (rho_r_out(mspin))
874  CALL qs_rho_set(rho_output, rho_r=rho_r_out)
875  IF (mspin > nspins) THEN
876  DO i = 1, mspin
877  CALL auxbas_pw_pool%create_pw(rho_r_out(i))
878  CALL pw_copy(rho_r_in(1), rho_r_out(i))
879  CALL pw_scale(rho_r_out(i), ospin)
880  END DO
881  ELSE
882  DO i = 1, nspins
883  CALL auxbas_pw_pool%create_pw(rho_r_out(i))
884  CALL pw_copy(rho_r_in(i), rho_r_out(i))
885  END DO
886  END IF
887  END IF
888 
889  ! rho_g
890  IF (ASSOCIATED(rho_g_in)) THEN
891  nspins = SIZE(rho_g_in)
892  cpassert(mspin >= nspins)
893  ALLOCATE (rho_g_out(mspin))
894  CALL qs_rho_set(rho_output, rho_g=rho_g_out)
895  IF (mspin > nspins) THEN
896  DO i = 1, mspin
897  CALL auxbas_pw_pool%create_pw(rho_g_out(i))
898  CALL pw_copy(rho_g_in(1), rho_g_out(i))
899  CALL pw_scale(rho_g_out(i), ospin)
900  END DO
901  ELSE
902  DO i = 1, nspins
903  CALL auxbas_pw_pool%create_pw(rho_g_out(i))
904  CALL pw_copy(rho_g_in(i), rho_g_out(i))
905  END DO
906  END IF
907  END IF
908 
909  ! SCCS
910  IF (ASSOCIATED(rho_r_sccs_in)) THEN
911  CALL qs_rho_set(rho_output, rho_r_sccs=rho_r_sccs_out)
912  CALL auxbas_pw_pool%create_pw(rho_r_sccs_out)
913  CALL pw_copy(rho_r_sccs_in, rho_r_sccs_out)
914  END IF
915 
916  ! drho_r
917  IF (ASSOCIATED(drho_r_in)) THEN
918  nspins = SIZE(drho_r_in)
919  cpassert(mspin >= nspins)
920  ALLOCATE (drho_r_out(3, mspin))
921  CALL qs_rho_set(rho_output, drho_r=drho_r_out)
922  IF (mspin > nspins) THEN
923  DO j = 1, mspin
924  DO i = 1, 3
925  CALL auxbas_pw_pool%create_pw(drho_r_out(i, j))
926  CALL pw_copy(drho_r_in(i, 1), drho_r_out(i, j))
927  CALL pw_scale(drho_r_out(i, j), ospin)
928  END DO
929  END DO
930  ELSE
931  DO j = 1, nspins
932  DO i = 1, 3
933  CALL auxbas_pw_pool%create_pw(drho_r_out(i, j))
934  CALL pw_copy(drho_r_in(i, j), drho_r_out(i, j))
935  END DO
936  END DO
937  END IF
938  END IF
939 
940  ! drho_g
941  IF (ASSOCIATED(drho_g_in)) THEN
942  nspins = SIZE(drho_g_in)
943  cpassert(mspin >= nspins)
944  ALLOCATE (drho_g_out(3, mspin))
945  CALL qs_rho_set(rho_output, drho_g=drho_g_out)
946  IF (mspin > nspins) THEN
947  DO j = 1, mspin
948  DO i = 1, 3
949  CALL auxbas_pw_pool%create_pw(drho_g_out(i, j))
950  CALL pw_copy(drho_g_in(i, 1), drho_g_out(i, j))
951  CALL pw_scale(drho_g_out(i, j), ospin)
952  END DO
953  END DO
954  ELSE
955  DO j = 1, nspins
956  DO i = 1, 3
957  CALL auxbas_pw_pool%create_pw(drho_g_out(i, j))
958  CALL pw_copy(drho_g_in(i, j), drho_g_out(i, j))
959  END DO
960  END DO
961  END IF
962  END IF
963 
964  ! tau_r
965  IF (ASSOCIATED(tau_r_in)) THEN
966  nspins = SIZE(tau_r_in)
967  cpassert(mspin >= nspins)
968  ALLOCATE (tau_r_out(mspin))
969  CALL qs_rho_set(rho_output, tau_r=tau_r_out)
970  IF (mspin > nspins) THEN
971  DO i = 1, mspin
972  CALL auxbas_pw_pool%create_pw(tau_r_out(i))
973  CALL pw_copy(tau_r_in(1), tau_r_out(i))
974  CALL pw_scale(tau_r_out(i), ospin)
975  END DO
976  ELSE
977  DO i = 1, nspins
978  CALL auxbas_pw_pool%create_pw(tau_r_out(i))
979  CALL pw_copy(tau_r_in(i), tau_r_out(i))
980  END DO
981  END IF
982  END IF
983 
984  ! tau_g
985  IF (ASSOCIATED(tau_g_in)) THEN
986  nspins = SIZE(tau_g_in)
987  cpassert(mspin >= nspins)
988  ALLOCATE (tau_g_out(mspin))
989  CALL qs_rho_set(rho_output, tau_g=tau_g_out)
990  IF (mspin > nspins) THEN
991  DO i = 1, mspin
992  CALL auxbas_pw_pool%create_pw(tau_g_out(i))
993  CALL pw_copy(tau_g_in(1), tau_g_out(i))
994  CALL pw_scale(tau_g_out(i), ospin)
995  END DO
996  ELSE
997  DO i = 1, nspins
998  CALL auxbas_pw_pool%create_pw(tau_g_out(i))
999  CALL pw_copy(tau_g_in(i), tau_g_out(i))
1000  END DO
1001  END IF
1002  END IF
1003 
1004  ! tot_rho_r
1005  IF (ASSOCIATED(tot_rho_r_in)) THEN
1006  nspins = SIZE(tot_rho_r_in)
1007  cpassert(mspin >= nspins)
1008  ALLOCATE (tot_rho_r_out(mspin))
1009  CALL qs_rho_set(rho_output, tot_rho_r=tot_rho_r_out)
1010  IF (mspin > nspins) THEN
1011  DO i = 1, mspin
1012  tot_rho_r_out(i) = tot_rho_r_in(1)*ospin
1013  END DO
1014  ELSE
1015  DO i = 1, nspins
1016  tot_rho_r_out(i) = tot_rho_r_in(i)
1017  END DO
1018  END IF
1019  END IF
1020 
1021  ! tot_rho_g
1022  IF (ASSOCIATED(tot_rho_g_in)) THEN
1023  nspins = SIZE(tot_rho_g_in)
1024  cpassert(mspin >= nspins)
1025  ALLOCATE (tot_rho_g_out(mspin))
1026  CALL qs_rho_set(rho_output, tot_rho_g=tot_rho_g_out)
1027  IF (mspin > nspins) THEN
1028  DO i = 1, mspin
1029  tot_rho_g_out(i) = tot_rho_g_in(1)*ospin
1030  END DO
1031  ELSE
1032  DO i = 1, nspins
1033  tot_rho_g_out(i) = tot_rho_g_in(i)
1034  END DO
1035  END IF
1036  END IF
1037 
1038  CALL qs_rho_set(rho_output, &
1039  rho_g_valid=rho_g_valid_in, &
1040  rho_r_valid=rho_r_valid_in, &
1041  drho_g_valid=drho_g_valid_in, &
1042  drho_r_valid=drho_r_valid_in, &
1043  tau_r_valid=tau_r_valid_in, &
1044  tau_g_valid=tau_g_valid_in, &
1045  soft_valid=soft_valid_in, &
1046  complex_rho_ao=complex_rho_ao)
1047 
1048  CALL timestop(handle)
1049 
1050  END SUBROUTINE qs_rho_copy
1051 
1052 ! **************************************************************************************************
1053 !> \brief rhoa = alpha*rhoa+beta*rhob
1054 !> \param rhoa ...
1055 !> \param rhob ...
1056 !> \param alpha ...
1057 !> \param beta ...
1058 ! **************************************************************************************************
1059  SUBROUTINE qs_rho_scale_and_add(rhoa, rhob, alpha, beta)
1060 
1061  TYPE(qs_rho_type), INTENT(IN) :: rhoa, rhob
1062  REAL(kind=dp), INTENT(IN) :: alpha, beta
1063 
1064  CHARACTER(len=*), PARAMETER :: routinen = 'qs_rho_scale_and_add'
1065 
1066  INTEGER :: handle, i, j, nspina, nspinb, nspins
1067  REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_g_a, tot_rho_g_b, tot_rho_r_a, &
1068  tot_rho_r_b
1069  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_a, rho_ao_b, rho_ao_im_a, &
1070  rho_ao_im_b
1071  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_a, rho_g_b, tau_g_a, tau_g_b
1072  TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER :: drho_g_a, drho_g_b
1073  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_a, rho_r_b, tau_r_a, tau_r_b
1074  TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER :: drho_r_a, drho_r_b
1075  TYPE(pw_r3d_rs_type), POINTER :: rho_r_sccs_a, rho_r_sccs_b
1076 
1077  CALL timeset(routinen, handle)
1078 
1079  NULLIFY (rho_ao_a, rho_ao_im_a, rho_r_a, rho_g_a, drho_r_a, &
1080  drho_g_a, tau_r_a, tau_g_a, tot_rho_r_a, tot_rho_g_a, rho_r_sccs_a)
1081 
1082  CALL qs_rho_get(rhoa, &
1083  rho_ao=rho_ao_a, &
1084  rho_ao_im=rho_ao_im_a, &
1085  rho_r=rho_r_a, &
1086  rho_g=rho_g_a, &
1087  drho_r=drho_r_a, &
1088  drho_g=drho_g_a, &
1089  tau_r=tau_r_a, &
1090  tau_g=tau_g_a, &
1091  tot_rho_r=tot_rho_r_a, &
1092  tot_rho_g=tot_rho_g_a, &
1093  rho_r_sccs=rho_r_sccs_a)
1094 
1095  NULLIFY (rho_ao_b, rho_ao_im_b, rho_r_b, rho_g_b, drho_r_b, &
1096  drho_g_b, tau_r_b, tau_g_b, tot_rho_r_b, tot_rho_g_b, rho_r_sccs_b)
1097 
1098  CALL qs_rho_get(rhob, &
1099  rho_ao=rho_ao_b, &
1100  rho_ao_im=rho_ao_im_b, &
1101  rho_r=rho_r_b, &
1102  rho_g=rho_g_b, &
1103  drho_r=drho_r_b, &
1104  drho_g=drho_g_b, &
1105  tau_r=tau_r_b, &
1106  tau_g=tau_g_b, &
1107  tot_rho_r=tot_rho_r_b, &
1108  tot_rho_g=tot_rho_g_b, &
1109  rho_r_sccs=rho_r_sccs_b)
1110  ! rho_ao
1111  IF (ASSOCIATED(rho_ao_a) .AND. ASSOCIATED(rho_ao_b)) THEN
1112  nspina = SIZE(rho_ao_a)
1113  nspinb = SIZE(rho_ao_b)
1114  nspins = min(nspina, nspinb)
1115  DO i = 1, nspins
1116  CALL dbcsr_add(rho_ao_a(i)%matrix, rho_ao_b(i)%matrix, alpha, beta)
1117  END DO
1118  END IF
1119 
1120  ! rho_ao_im
1121  IF (ASSOCIATED(rho_ao_im_a) .AND. ASSOCIATED(rho_ao_im_b)) THEN
1122  nspina = SIZE(rho_ao_im_a)
1123  nspinb = SIZE(rho_ao_im_b)
1124  nspins = min(nspina, nspinb)
1125  DO i = 1, nspins
1126  CALL dbcsr_add(rho_ao_im_a(i)%matrix, rho_ao_im_b(i)%matrix, alpha, beta)
1127  END DO
1128  END IF
1129 
1130  ! rho_r
1131  IF (ASSOCIATED(rho_r_a) .AND. ASSOCIATED(rho_r_b)) THEN
1132  nspina = SIZE(rho_ao_a)
1133  nspinb = SIZE(rho_ao_b)
1134  nspins = min(nspina, nspinb)
1135  DO i = 1, nspins
1136  CALL pw_axpy(rho_r_b(i), rho_r_a(i), beta, alpha)
1137  END DO
1138  END IF
1139 
1140  ! rho_g
1141  IF (ASSOCIATED(rho_g_a) .AND. ASSOCIATED(rho_g_b)) THEN
1142  nspina = SIZE(rho_ao_a)
1143  nspinb = SIZE(rho_ao_b)
1144  nspins = min(nspina, nspinb)
1145  DO i = 1, nspins
1146  CALL pw_axpy(rho_g_b(i), rho_g_a(i), beta, alpha)
1147  END DO
1148  END IF
1149 
1150  ! SCCS
1151  IF (ASSOCIATED(rho_r_sccs_a) .AND. ASSOCIATED(rho_r_sccs_b)) THEN
1152  CALL pw_axpy(rho_r_sccs_b, rho_r_sccs_a, beta, alpha)
1153  END IF
1154 
1155  ! drho_r
1156  IF (ASSOCIATED(drho_r_a) .AND. ASSOCIATED(drho_r_b)) THEN
1157  cpassert(all(shape(drho_r_a) == shape(drho_r_b))) ! not implemented
1158  DO j = 1, SIZE(drho_r_a, 2)
1159  DO i = 1, SIZE(drho_r_a, 1)
1160  CALL pw_axpy(drho_r_b(i, j), drho_r_a(i, j), beta, alpha)
1161  END DO
1162  END DO
1163  END IF
1164 
1165  ! drho_g
1166  IF (ASSOCIATED(drho_g_a) .AND. ASSOCIATED(drho_g_b)) THEN
1167  cpassert(all(shape(drho_g_a) == shape(drho_g_b))) ! not implemented
1168  DO j = 1, SIZE(drho_g_a, 2)
1169  DO i = 1, SIZE(drho_g_a, 1)
1170  CALL pw_axpy(drho_g_b(i, j), drho_g_a(i, j), beta, alpha)
1171  END DO
1172  END DO
1173  END IF
1174 
1175  ! tau_r
1176  IF (ASSOCIATED(tau_r_a) .AND. ASSOCIATED(tau_r_b)) THEN
1177  nspina = SIZE(rho_ao_a)
1178  nspinb = SIZE(rho_ao_b)
1179  nspins = min(nspina, nspinb)
1180  DO i = 1, nspins
1181  CALL pw_axpy(tau_r_b(i), tau_r_a(i), beta, alpha)
1182  END DO
1183  END IF
1184 
1185  ! tau_g
1186  IF (ASSOCIATED(tau_g_a) .AND. ASSOCIATED(tau_g_b)) THEN
1187  nspina = SIZE(rho_ao_a)
1188  nspinb = SIZE(rho_ao_b)
1189  nspins = min(nspina, nspinb)
1190  DO i = 1, nspins
1191  CALL pw_axpy(tau_g_b(i), tau_g_a(i), beta, alpha)
1192  END DO
1193  END IF
1194 
1195  ! tot_rho_r
1196  IF (ASSOCIATED(tot_rho_r_a) .AND. ASSOCIATED(tot_rho_r_b)) THEN
1197  nspina = SIZE(rho_ao_a)
1198  nspinb = SIZE(rho_ao_b)
1199  nspins = min(nspina, nspinb)
1200  DO i = 1, nspins
1201  tot_rho_r_a(i) = alpha*tot_rho_r_a(i) + beta*tot_rho_r_b(i)
1202  END DO
1203  END IF
1204 
1205  ! tot_rho_g
1206  IF (ASSOCIATED(tot_rho_g_a) .AND. ASSOCIATED(tot_rho_g_b)) THEN
1207  nspina = SIZE(rho_ao_a)
1208  nspinb = SIZE(rho_ao_b)
1209  nspins = min(nspina, nspinb)
1210  DO i = 1, nspins
1211  tot_rho_g_a(i) = alpha*tot_rho_g_a(i) + beta*tot_rho_g_b(i)
1212  END DO
1213  END IF
1214 
1215  CALL timestop(handle)
1216 
1217  END SUBROUTINE qs_rho_scale_and_add
1218 
1219 ! **************************************************************************************************
1220 !> \brief Duplicates a pointer physically
1221 !> \param rho_input The rho structure to be duplicated
1222 !> \param rho_output The duplicate rho structure
1223 !> \param qs_env The QS environment from which the auxiliary PW basis-set
1224 !> pool is taken
1225 !> \par History
1226 !> 07.2005 initial create [tdk]
1227 !> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch)
1228 !> \note
1229 !> Associated pointers are deallocated, nullified pointers are NOT accepted!
1230 ! **************************************************************************************************
1231  SUBROUTINE duplicate_rho_type(rho_input, rho_output, qs_env)
1232 
1233  TYPE(qs_rho_type), INTENT(INOUT) :: rho_input, rho_output
1234  TYPE(qs_environment_type), POINTER :: qs_env
1235 
1236  CHARACTER(len=*), PARAMETER :: routinen = 'duplicate_rho_type'
1237 
1238  INTEGER :: handle, i, j, nspins
1239  LOGICAL :: complex_rho_ao_in, drho_g_valid_in, drho_r_valid_in, rho_g_valid_in, &
1240  rho_r_valid_in, soft_valid_in, tau_g_valid_in, tau_r_valid_in
1241  REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_g_in, tot_rho_g_out, &
1242  tot_rho_r_in, tot_rho_r_out
1243  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_im_in, rho_ao_im_out, rho_ao_in, &
1244  rho_ao_out
1245  TYPE(dft_control_type), POINTER :: dft_control
1246  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_in, rho_g_out, tau_g_in, tau_g_out
1247  TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER :: drho_g_in, drho_g_out
1248  TYPE(pw_env_type), POINTER :: pw_env
1249  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1250  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_in, rho_r_out, tau_r_in, tau_r_out
1251  TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER :: drho_r_in, drho_r_out
1252  TYPE(pw_r3d_rs_type), POINTER :: rho_r_sccs_in, rho_r_sccs_out
1253 
1254  CALL timeset(routinen, handle)
1255 
1256  NULLIFY (dft_control, pw_env, auxbas_pw_pool)
1257  NULLIFY (rho_ao_in, rho_ao_out, rho_ao_im_in, rho_ao_im_out)
1258  NULLIFY (rho_r_in, rho_r_out, rho_g_in, rho_g_out, drho_r_in, drho_r_out)
1259  NULLIFY (drho_g_in, drho_g_out, tau_r_in, tau_r_out, tau_g_in, tau_g_out)
1260  NULLIFY (tot_rho_r_in, tot_rho_r_out, tot_rho_g_in, tot_rho_g_out)
1261  NULLIFY (rho_r_sccs_in, rho_r_sccs_out)
1262 
1263  cpassert(ASSOCIATED(qs_env))
1264 
1265  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, dft_control=dft_control)
1266  CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool)
1267  nspins = dft_control%nspins
1268 
1269  CALL qs_rho_clear(rho_output)
1270 
1271  CALL qs_rho_get(rho_input, &
1272  rho_ao=rho_ao_in, &
1273  rho_ao_im=rho_ao_im_in, &
1274  rho_r=rho_r_in, &
1275  rho_g=rho_g_in, &
1276  drho_r=drho_r_in, &
1277  drho_g=drho_g_in, &
1278  tau_r=tau_r_in, &
1279  tau_g=tau_g_in, &
1280  tot_rho_r=tot_rho_r_in, &
1281  tot_rho_g=tot_rho_g_in, &
1282  rho_g_valid=rho_g_valid_in, &
1283  rho_r_valid=rho_r_valid_in, &
1284  drho_g_valid=drho_g_valid_in, &
1285  drho_r_valid=drho_r_valid_in, &
1286  tau_r_valid=tau_r_valid_in, &
1287  tau_g_valid=tau_g_valid_in, &
1288  rho_r_sccs=rho_r_sccs_in, &
1289  soft_valid=soft_valid_in, &
1290  complex_rho_ao=complex_rho_ao_in)
1291 
1292  ! rho_ao
1293  IF (ASSOCIATED(rho_ao_in)) THEN
1294  CALL dbcsr_allocate_matrix_set(rho_ao_out, nspins)
1295  CALL qs_rho_set(rho_output, rho_ao=rho_ao_out)
1296  DO i = 1, nspins
1297  ALLOCATE (rho_ao_out(i)%matrix)
1298  CALL dbcsr_copy(rho_ao_out(i)%matrix, rho_ao_in(i)%matrix, &
1299  name="myDensityMatrix_for_Spin_"//trim(adjustl(cp_to_string(i))))
1300  CALL dbcsr_set(rho_ao_out(i)%matrix, 0.0_dp)
1301  END DO
1302  END IF
1303 
1304  ! rho_ao_im
1305  IF (ASSOCIATED(rho_ao_im_in)) THEN
1306  CALL dbcsr_allocate_matrix_set(rho_ao_im_out, nspins)
1307  CALL qs_rho_set(rho_output, rho_ao=rho_ao_im_out)
1308  DO i = 1, nspins
1309  ALLOCATE (rho_ao_im_out(i)%matrix)
1310  CALL dbcsr_copy(rho_ao_im_out(i)%matrix, rho_ao_im_in(i)%matrix, &
1311  name="myImagDensityMatrix_for_Spin_"//trim(adjustl(cp_to_string(i))))
1312  CALL dbcsr_set(rho_ao_im_out(i)%matrix, 0.0_dp)
1313  END DO
1314  END IF
1315 
1316  ! rho_r
1317  IF (ASSOCIATED(rho_r_in)) THEN
1318  ALLOCATE (rho_r_out(nspins))
1319  CALL qs_rho_set(rho_output, rho_r=rho_r_out)
1320  DO i = 1, nspins
1321  CALL auxbas_pw_pool%create_pw(rho_r_out(i))
1322  CALL pw_copy(rho_r_in(i), rho_r_out(i))
1323  END DO
1324  END IF
1325 
1326  ! rho_g
1327  IF (ASSOCIATED(rho_g_in)) THEN
1328  ALLOCATE (rho_g_out(nspins))
1329  CALL qs_rho_set(rho_output, rho_g=rho_g_out)
1330  DO i = 1, nspins
1331  CALL auxbas_pw_pool%create_pw(rho_g_out(i))
1332  CALL pw_copy(rho_g_in(i), rho_g_out(i))
1333  END DO
1334  END IF
1335 
1336  ! SCCS
1337  IF (ASSOCIATED(rho_r_sccs_in)) THEN
1338  CALL qs_rho_set(rho_output, rho_r_sccs=rho_r_sccs_out)
1339  CALL auxbas_pw_pool%create_pw(rho_r_sccs_out)
1340  CALL pw_copy(rho_r_sccs_in, rho_r_sccs_out)
1341  END IF
1342 
1343  ! drho_r and drho_g are only needed if calculated by collocation
1344  IF (dft_control%drho_by_collocation) THEN
1345  ! drho_r
1346  IF (ASSOCIATED(drho_r_in)) THEN
1347  ALLOCATE (drho_r_out(3, nspins))
1348  CALL qs_rho_set(rho_output, drho_r=drho_r_out)
1349  DO j = 1, nspins
1350  DO i = 1, 3
1351  CALL auxbas_pw_pool%create_pw(drho_r_out(i, j))
1352  CALL pw_copy(drho_r_in(i, j), drho_r_out(i, j))
1353  END DO
1354  END DO
1355  END IF
1356 
1357  ! drho_g
1358  IF (ASSOCIATED(drho_g_in)) THEN
1359  ALLOCATE (drho_g_out(3, nspins))
1360  CALL qs_rho_set(rho_output, drho_g=drho_g_out)
1361  DO j = 1, nspins
1362  DO i = 1, 3
1363  CALL auxbas_pw_pool%create_pw(drho_g_out(i, j))
1364  CALL pw_copy(drho_g_in(i, j), drho_g_out(i, j))
1365  END DO
1366  END DO
1367  END IF
1368  END IF
1369 
1370  ! tau_r and tau_g are only needed in the case of Meta-GGA XC-functionals
1371  ! are used. Therefore they are only allocated if
1372  ! dft_control%use_kinetic_energy_density is true
1373  IF (dft_control%use_kinetic_energy_density) THEN
1374  ! tau_r
1375  IF (ASSOCIATED(tau_r_in)) THEN
1376  ALLOCATE (tau_r_out(nspins))
1377  CALL qs_rho_set(rho_output, tau_r=tau_r_out)
1378  DO i = 1, nspins
1379  CALL auxbas_pw_pool%create_pw(tau_r_out(i))
1380  CALL pw_copy(tau_r_in(i), tau_r_out(i))
1381  END DO
1382  END IF
1383 
1384  ! tau_g
1385  IF (ASSOCIATED(tau_g_in)) THEN
1386  ALLOCATE (tau_g_out(nspins))
1387  CALL qs_rho_set(rho_output, tau_g=tau_g_out)
1388  DO i = 1, nspins
1389  CALL auxbas_pw_pool%create_pw(tau_g_out(i))
1390  CALL pw_copy(tau_g_in(i), tau_g_out(i))
1391  END DO
1392  END IF
1393  END IF
1394 
1395  CALL qs_rho_set(rho_output, &
1396  rho_g_valid=rho_g_valid_in, &
1397  rho_r_valid=rho_r_valid_in, &
1398  drho_g_valid=drho_g_valid_in, &
1399  drho_r_valid=drho_r_valid_in, &
1400  tau_r_valid=tau_r_valid_in, &
1401  tau_g_valid=tau_g_valid_in, &
1402  soft_valid=soft_valid_in, &
1403  complex_rho_ao=complex_rho_ao_in)
1404 
1405  ! tot_rho_r
1406  IF (ASSOCIATED(tot_rho_r_in)) THEN
1407  ALLOCATE (tot_rho_r_out(nspins))
1408  CALL qs_rho_set(rho_output, tot_rho_r=tot_rho_r_out)
1409  DO i = 1, nspins
1410  tot_rho_r_out(i) = tot_rho_r_in(i)
1411  END DO
1412  END IF
1413 
1414  ! tot_rho_g
1415  IF (ASSOCIATED(tot_rho_g_in)) THEN
1416  ALLOCATE (tot_rho_g_out(nspins))
1417  CALL qs_rho_set(rho_output, tot_rho_g=tot_rho_g_out)
1418  DO i = 1, nspins
1419  tot_rho_g_out(i) = tot_rho_g_in(i)
1420  END DO
1421 
1422  END IF
1423 
1424  CALL timestop(handle)
1425 
1426  END SUBROUTINE duplicate_rho_type
1427 
1428 ! **************************************************************************************************
1429 !> \brief (Re-)allocates rho_ao_im from real part rho_ao
1430 !> \param rho ...
1431 !> \param qs_env ...
1432 ! **************************************************************************************************
1433  SUBROUTINE allocate_rho_ao_imag_from_real(rho, qs_env)
1434  TYPE(qs_rho_type), POINTER :: rho
1435  TYPE(qs_environment_type), POINTER :: qs_env
1436 
1437  CHARACTER(LEN=default_string_length) :: headline
1438  INTEGER :: i, ic, nimages, nspins
1439  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_im_kp, rho_ao_kp
1440  TYPE(dbcsr_type), POINTER :: template
1441  TYPE(dft_control_type), POINTER :: dft_control
1442  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1443  POINTER :: sab_orb
1444 
1445  NULLIFY (rho_ao_im_kp, rho_ao_kp, dft_control, template, sab_orb)
1446 
1447  CALL get_qs_env(qs_env, &
1448  dft_control=dft_control, &
1449  sab_orb=sab_orb)
1450 
1451  CALL qs_rho_get(rho, rho_ao_im_kp=rho_ao_im_kp, rho_ao_kp=rho_ao_kp)
1452 
1453  nspins = dft_control%nspins
1454  nimages = dft_control%nimages
1455 
1456  cpassert(nspins .EQ. SIZE(rho_ao_kp, 1))
1457  cpassert(nimages .EQ. SIZE(rho_ao_kp, 2))
1458 
1459  CALL dbcsr_allocate_matrix_set(rho_ao_im_kp, nspins, nimages)
1460  CALL qs_rho_set(rho, rho_ao_im_kp=rho_ao_im_kp)
1461  DO i = 1, nspins
1462  DO ic = 1, nimages
1463  IF (nspins > 1) THEN
1464  IF (i == 1) THEN
1465  headline = "IMAGINARY PART OF DENSITY MATRIX FOR ALPHA SPIN"
1466  ELSE
1467  headline = "IMAGINARY PART OF DENSITY MATRIX FOR BETA SPIN"
1468  END IF
1469  ELSE
1470  headline = "IMAGINARY PART OF DENSITY MATRIX"
1471  END IF
1472  ALLOCATE (rho_ao_im_kp(i, ic)%matrix)
1473  template => rho_ao_kp(i, ic)%matrix ! base on real part, but anti-symmetric
1474  CALL dbcsr_create(matrix=rho_ao_im_kp(i, ic)%matrix, template=template, &
1475  name=trim(headline), matrix_type=dbcsr_type_antisymmetric, nze=0)
1476  CALL cp_dbcsr_alloc_block_from_nbl(rho_ao_im_kp(i, ic)%matrix, sab_orb)
1477  CALL dbcsr_set(rho_ao_im_kp(i, ic)%matrix, 0.0_dp)
1478  END DO
1479  END DO
1480 
1481  END SUBROUTINE allocate_rho_ao_imag_from_real
1482 
1483 END MODULE qs_rho_methods
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
admm_types
Types and set/get functions for auxiliary density matrix methods.
Definition: admm_types.F:15
admm_types::get_admm_env
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Definition: admm_types.F:593
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition: cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition: cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition: kpoint_types.F:333
lri_environment_methods
Calculates integral matrices for LRIGPW method lri : local resolution of the identity.
Definition: lri_environment_methods.F:18
lri_environment_methods::calculate_lri_densities
subroutine, public calculate_lri_densities(lri_env, lri_density, qs_env, pmatrix, cell_to_index, lri_rho_struct, atomic_kind_set, para_env, response_density)
performs the fitting of the density and distributes the fitted density on the grid
Definition: lri_environment_methods.F:569
lri_environment_types
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Definition: lri_environment_types.F:18
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
pw_env_types
container for various plainwaves related things
Definition: pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
pw_methods
Definition: pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
pw_types
Definition: pw_types.F:24
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition: qs_collocate_density.F:38
qs_collocate_density::calculate_drho_elec
subroutine, public calculate_drho_elec(matrix_p, matrix_p_kp, drho, drho_gspace, qs_env, soft_valid, basis_type)
computes the gradient of the density corresponding to a given density matrix on the grid
Definition: qs_collocate_density.F:1876
qs_collocate_density::calculate_rho_elec
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
Definition: qs_collocate_density.F:1710
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Set the QUICKSTEP environment.
Definition: qs_environment_types.F:1035
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_ks_types
Definition: qs_ks_types.F:15
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition: qs_ks_types.F:330
qs_local_rho_types
Definition: qs_local_rho_types.F:8
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_oce_types
Definition: qs_oce_types.F:8
qs_rho_atom_methods
Definition: qs_rho_atom_methods.F:7
qs_rho_atom_methods::calculate_rho_atom_coeff
subroutine, public calculate_rho_atom_coeff(qs_env, rho_ao, rho_atom_set, qs_kind_set, oce, sab, para_env)
...
Definition: qs_rho_atom_methods.F:420
qs_rho_atom_types
Definition: qs_rho_atom_types.F:8
qs_rho_methods
methods of the rho structure (defined in qs_rho_types)
Definition: qs_rho_methods.F:15
qs_rho_methods::qs_rho_update_tddfpt
subroutine, public qs_rho_update_tddfpt(rho_struct, qs_env, pw_env_external, task_list_external, para_env_external, tddfpt_lri_env, tddfpt_lri_density)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g
Definition: qs_rho_methods.F:670
qs_rho_methods::allocate_rho_ao_imag_from_real
subroutine, public allocate_rho_ao_imag_from_real(rho, qs_env)
(Re-)allocates rho_ao_im from real part rho_ao
Definition: qs_rho_methods.F:1434
qs_rho_methods::qs_rho_update_rho
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
Definition: qs_rho_methods.F:375
qs_rho_methods::duplicate_rho_type
subroutine, public duplicate_rho_type(rho_input, rho_output, qs_env)
Duplicates a pointer physically.
Definition: qs_rho_methods.F:1232
qs_rho_methods::qs_rho_copy
subroutine, public qs_rho_copy(rho_input, rho_output, auxbas_pw_pool, mspin)
Allocate a density structure and fill it with data from an input structure SIZE(rho_input) == mspin =...
Definition: qs_rho_methods.F:766
qs_rho_methods::qs_rho_scale_and_add
subroutine, public qs_rho_scale_and_add(rhoa, rhob, alpha, beta)
rhoa = alpha*rhoa+beta*rhob
Definition: qs_rho_methods.F:1060
qs_rho_methods::qs_rho_rebuild
subroutine, public qs_rho_rebuild(rho, qs_env, rebuild_ao, rebuild_grids, admm, pw_env_external)
rebuilds rho (if necessary allocating and initializing it)
Definition: qs_rho_methods.F:96
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_set
subroutine, public qs_rho_set(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
...
Definition: qs_rho_types.F:308
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
qs_rho_types::qs_rho_clear
subroutine, public qs_rho_clear(rho_struct)
Deallocates all components, without deallocating rho_struct itself.
Definition: qs_rho_types.F:125
ri_environment_methods
Calculates integral matrices for RIGPW method.
Definition: ri_environment_methods.F:18
ri_environment_methods::calculate_ri_densities
subroutine, public calculate_ri_densities(lri_env, qs_env, pmatrix, lri_rho_struct, atomic_kind_set, para_env)
performs the fitting of the density and distributes the fitted density on the grid
Definition: ri_environment_methods.F:358
task_list_types
types for task lists
Definition: task_list_types.F:14