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