Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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_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_transfer,&
40 : pw_zero
41 : USE pw_pool_types, ONLY: pw_pool_type
42 : USE pw_types, ONLY: pw_c1d_gs_type,&
43 : pw_r3d_rs_type
44 : USE qs_collocate_density, ONLY: calculate_drho_elec,&
45 : calculate_rho_elec
46 : USE qs_environment_types, ONLY: get_qs_env,&
47 : qs_environment_type,&
48 : set_qs_env
49 : USE qs_harris_types, ONLY: harris_type
50 : USE qs_harris_utils, ONLY: calculate_harris_density
51 : USE qs_kind_types, ONLY: qs_kind_type
52 : USE qs_ks_types, ONLY: get_ks_env,&
53 : qs_ks_env_type
54 : USE qs_local_rho_types, ONLY: local_rho_type
55 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
56 : USE qs_oce_types, ONLY: oce_matrix_type
57 : USE qs_rho_atom_methods, ONLY: calculate_rho_atom_coeff
58 : USE qs_rho_atom_types, ONLY: rho_atom_type
59 : USE qs_rho_types, ONLY: qs_rho_clear,&
60 : qs_rho_get,&
61 : qs_rho_set,&
62 : qs_rho_type
63 : USE ri_environment_methods, ONLY: calculate_ri_densities
64 : USE task_list_types, ONLY: task_list_type
65 : #include "./base/base_uses.f90"
66 :
67 : IMPLICIT NONE
68 : PRIVATE
69 :
70 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
71 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho_methods'
72 :
73 : PUBLIC :: qs_rho_update_rho, qs_rho_update_tddfpt, &
74 : qs_rho_rebuild, qs_rho_copy, qs_rho_scale_and_add, &
75 : qs_rho_scale_and_add_b, qs_rho_transfer
76 : PUBLIC :: duplicate_rho_type, allocate_rho_ao_imag_from_real
77 :
78 : CONTAINS
79 :
80 : ! **************************************************************************************************
81 : !> \brief rebuilds rho (if necessary allocating and initializing it)
82 : !> \param rho the rho type to rebuild (defaults to qs_env%rho)
83 : !> \param qs_env the environment to which rho belongs
84 : !> \param rebuild_ao if it is necessary to rebuild rho_ao. Defaults to true.
85 : !> \param rebuild_grids if it in necessary to rebuild rho_r and rho_g.
86 : !> Defaults to false.
87 : !> \param admm (use aux_fit basis)
88 : !> \param pw_env_external external plane wave environment
89 : !> \par History
90 : !> 11.2002 created replacing qs_rho_create and qs_env_rebuild_rho[fawzi]
91 : !> \author Fawzi Mohamed
92 : !> \note
93 : !> needs updated pw pools, s, s_mstruct and h in qs_env.
94 : !> The use of p to keep the structure of h (needed for the forces)
95 : !> is ugly and should be removed.
96 : !> Change so that it does not allocate a subcomponent if it is not
97 : !> associated and not requested?
98 : ! **************************************************************************************************
99 71108 : SUBROUTINE qs_rho_rebuild(rho, qs_env, rebuild_ao, rebuild_grids, admm, pw_env_external)
100 : TYPE(qs_rho_type), INTENT(INOUT) :: rho
101 : TYPE(qs_environment_type), POINTER :: qs_env
102 : LOGICAL, INTENT(in), OPTIONAL :: rebuild_ao, rebuild_grids, admm
103 : TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
104 :
105 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_rho_rebuild'
106 :
107 : CHARACTER(LEN=default_string_length) :: headline
108 : INTEGER :: handle, i, ic, j, nimg, nspins
109 : LOGICAL :: do_kpoints, my_admm, my_rebuild_ao, &
110 : my_rebuild_grids, rho_ao_is_complex
111 35554 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
112 35554 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, rho_ao_im_kp, rho_ao_kp
113 : TYPE(dbcsr_type), POINTER :: refmatrix, tmatrix
114 : TYPE(dft_control_type), POINTER :: dft_control
115 : TYPE(kpoint_type), POINTER :: kpoints
116 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
117 35554 : POINTER :: sab_orb
118 35554 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, tau_g
119 35554 : TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER :: drho_g
120 : TYPE(pw_env_type), POINTER :: pw_env
121 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
122 35554 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau_r
123 35554 : TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER :: drho_r
124 : TYPE(pw_r3d_rs_type), POINTER :: rho_r_sccs
125 :
126 35554 : CALL timeset(routineN, handle)
127 :
128 35554 : NULLIFY (pw_env, auxbas_pw_pool, matrix_s_kp, dft_control)
129 35554 : NULLIFY (tot_rho_r, rho_ao_kp, rho_r, rho_g, drho_r, drho_g, tau_r, tau_g, rho_ao_im_kp)
130 35554 : NULLIFY (rho_r_sccs)
131 35554 : NULLIFY (sab_orb)
132 35554 : my_rebuild_ao = .TRUE.
133 35554 : my_rebuild_grids = .TRUE.
134 35554 : my_admm = .FALSE.
135 35554 : IF (PRESENT(rebuild_ao)) my_rebuild_ao = rebuild_ao
136 35554 : IF (PRESENT(rebuild_grids)) my_rebuild_grids = rebuild_grids
137 35554 : IF (PRESENT(admm)) my_admm = admm
138 :
139 : CALL get_qs_env(qs_env, &
140 : kpoints=kpoints, &
141 : do_kpoints=do_kpoints, &
142 : pw_env=pw_env, &
143 35554 : dft_control=dft_control)
144 35554 : IF (PRESENT(pw_env_external)) THEN
145 1108 : pw_env => pw_env_external
146 : END IF
147 :
148 35554 : nimg = dft_control%nimages
149 :
150 35554 : IF (my_admm) THEN
151 2172 : CALL get_admm_env(qs_env%admm_env, sab_aux_fit=sab_orb, matrix_s_aux_fit_kp=matrix_s_kp)
152 : ELSE
153 33382 : CALL get_qs_env(qs_env, matrix_s_kp=matrix_s_kp)
154 :
155 33382 : IF (do_kpoints) THEN
156 3826 : CALL get_kpoint_info(kpoints, sab_nl=sab_orb)
157 : ELSE
158 29556 : CALL get_qs_env(qs_env, sab_orb=sab_orb)
159 : END IF
160 : END IF
161 35554 : refmatrix => matrix_s_kp(1, 1)%matrix
162 :
163 35554 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
164 35554 : nspins = dft_control%nspins
165 :
166 : CALL qs_rho_get(rho, &
167 : tot_rho_r=tot_rho_r, &
168 : rho_ao_kp=rho_ao_kp, &
169 : rho_ao_im_kp=rho_ao_im_kp, &
170 : rho_r=rho_r, &
171 : rho_g=rho_g, &
172 : drho_r=drho_r, &
173 : drho_g=drho_g, &
174 : tau_r=tau_r, &
175 : tau_g=tau_g, &
176 : rho_r_sccs=rho_r_sccs, &
177 35554 : complex_rho_ao=rho_ao_is_complex)
178 :
179 35554 : IF (.NOT. ASSOCIATED(tot_rho_r)) THEN
180 43476 : ALLOCATE (tot_rho_r(nspins))
181 31679 : tot_rho_r = 0.0_dp
182 14492 : CALL qs_rho_set(rho, tot_rho_r=tot_rho_r)
183 : END IF
184 :
185 : ! rho_ao
186 35554 : IF (my_rebuild_ao .OR. (.NOT. ASSOCIATED(rho_ao_kp))) THEN
187 33582 : IF (ASSOCIATED(rho_ao_kp)) THEN
188 21028 : CALL dbcsr_deallocate_matrix_set(rho_ao_kp)
189 : END IF
190 : ! Create a new density matrix set
191 33582 : CALL dbcsr_allocate_matrix_set(rho_ao_kp, nspins, nimg)
192 33582 : CALL qs_rho_set(rho, rho_ao_kp=rho_ao_kp)
193 71854 : DO i = 1, nspins
194 423082 : DO ic = 1, nimg
195 351228 : IF (nspins > 1) THEN
196 77012 : IF (i == 1) THEN
197 38506 : headline = "DENSITY MATRIX FOR ALPHA SPIN"
198 : ELSE
199 38506 : headline = "DENSITY MATRIX FOR BETA SPIN"
200 : END IF
201 : ELSE
202 274216 : headline = "DENSITY MATRIX"
203 : END IF
204 351228 : ALLOCATE (rho_ao_kp(i, ic)%matrix)
205 351228 : tmatrix => rho_ao_kp(i, ic)%matrix
206 : CALL dbcsr_create(matrix=tmatrix, template=refmatrix, name=TRIM(headline), &
207 351228 : matrix_type=dbcsr_type_symmetric)
208 351228 : CALL cp_dbcsr_alloc_block_from_nbl(tmatrix, sab_orb)
209 389500 : CALL dbcsr_set(tmatrix, 0.0_dp)
210 : END DO
211 : END DO
212 33582 : IF (rho_ao_is_complex) THEN
213 378 : IF (ASSOCIATED(rho_ao_im_kp)) THEN
214 378 : CALL dbcsr_deallocate_matrix_set(rho_ao_im_kp)
215 : END IF
216 378 : CALL dbcsr_allocate_matrix_set(rho_ao_im_kp, nspins, nimg)
217 378 : CALL qs_rho_set(rho, rho_ao_im_kp=rho_ao_im_kp)
218 836 : DO i = 1, nspins
219 1294 : DO ic = 1, nimg
220 458 : IF (nspins > 1) THEN
221 160 : IF (i == 1) THEN
222 80 : headline = "IMAGINARY PART OF DENSITY MATRIX FOR ALPHA SPIN"
223 : ELSE
224 80 : headline = "IMAGINARY PART OF DENSITY MATRIX FOR BETA SPIN"
225 : END IF
226 : ELSE
227 298 : headline = "IMAGINARY PART OF DENSITY MATRIX"
228 : END IF
229 458 : ALLOCATE (rho_ao_im_kp(i, ic)%matrix)
230 458 : tmatrix => rho_ao_im_kp(i, ic)%matrix
231 : CALL dbcsr_create(matrix=tmatrix, template=refmatrix, name=TRIM(headline), &
232 458 : matrix_type=dbcsr_type_antisymmetric)
233 458 : CALL cp_dbcsr_alloc_block_from_nbl(tmatrix, sab_orb)
234 916 : CALL dbcsr_set(tmatrix, 0.0_dp)
235 : END DO
236 : END DO
237 : END IF
238 : END IF
239 :
240 : ! rho_r
241 35554 : IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(rho_r)) THEN
242 35554 : IF (ASSOCIATED(rho_r)) THEN
243 44259 : DO i = 1, SIZE(rho_r)
244 44259 : CALL rho_r(i)%release()
245 : END DO
246 21062 : DEALLOCATE (rho_r)
247 : END IF
248 147046 : ALLOCATE (rho_r(nspins))
249 35554 : CALL qs_rho_set(rho, rho_r=rho_r)
250 75938 : DO i = 1, nspins
251 75938 : CALL auxbas_pw_pool%create_pw(rho_r(i))
252 : END DO
253 : END IF
254 :
255 : ! rho_g
256 35554 : IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(rho_g)) THEN
257 35554 : IF (ASSOCIATED(rho_g)) THEN
258 44259 : DO i = 1, SIZE(rho_g)
259 44259 : CALL rho_g(i)%release()
260 : END DO
261 21062 : DEALLOCATE (rho_g)
262 : END IF
263 147046 : ALLOCATE (rho_g(nspins))
264 35554 : CALL qs_rho_set(rho, rho_g=rho_g)
265 75938 : DO i = 1, nspins
266 75938 : CALL auxbas_pw_pool%create_pw(rho_g(i))
267 : END DO
268 : END IF
269 :
270 : ! SCCS
271 35554 : IF (dft_control%do_sccs) THEN
272 14 : IF (my_rebuild_grids .OR. (.NOT. ASSOCIATED(rho_r_sccs))) THEN
273 14 : IF (ASSOCIATED(rho_r_sccs)) THEN
274 2 : CALL rho_r_sccs%release()
275 2 : DEALLOCATE (rho_r_sccs)
276 : END IF
277 14 : ALLOCATE (rho_r_sccs)
278 14 : CALL qs_rho_set(rho, rho_r_sccs=rho_r_sccs)
279 14 : CALL auxbas_pw_pool%create_pw(rho_r_sccs)
280 14 : CALL pw_zero(rho_r_sccs)
281 : END IF
282 : END IF
283 :
284 : ! allocate drho_r and drho_g if xc_deriv_collocate
285 35554 : IF (dft_control%drho_by_collocation) THEN
286 : ! drho_r
287 0 : IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(drho_r)) THEN
288 0 : IF (ASSOCIATED(drho_r)) THEN
289 0 : DO j = 1, SIZE(drho_r, 2)
290 0 : DO i = 1, SIZE(drho_r, 1)
291 0 : CALL drho_r(i, j)%release()
292 : END DO
293 : END DO
294 0 : DEALLOCATE (drho_r)
295 : END IF
296 0 : ALLOCATE (drho_r(3, nspins))
297 0 : CALL qs_rho_set(rho, drho_r=drho_r)
298 0 : DO j = 1, nspins
299 0 : DO i = 1, 3
300 0 : CALL auxbas_pw_pool%create_pw(drho_r(i, j))
301 : END DO
302 : END DO
303 : END IF
304 : ! drho_g
305 0 : IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(drho_g)) THEN
306 0 : IF (ASSOCIATED(drho_g)) THEN
307 0 : DO j = 1, SIZE(drho_g, 2)
308 0 : DO i = 1, SIZE(drho_r, 1)
309 0 : CALL drho_g(i, j)%release()
310 : END DO
311 : END DO
312 0 : DEALLOCATE (drho_g)
313 : END IF
314 0 : ALLOCATE (drho_g(3, nspins))
315 0 : CALL qs_rho_set(rho, drho_g=drho_g)
316 0 : DO j = 1, nspins
317 0 : DO i = 1, 3
318 0 : CALL auxbas_pw_pool%create_pw(drho_g(i, j))
319 : END DO
320 : END DO
321 : END IF
322 : END IF
323 :
324 : ! allocate tau_r and tau_g if use_kinetic_energy_density
325 35554 : IF (dft_control%use_kinetic_energy_density) THEN
326 : ! tau_r
327 780 : IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(tau_r)) THEN
328 780 : IF (ASSOCIATED(tau_r)) THEN
329 458 : DO i = 1, SIZE(tau_r)
330 458 : CALL tau_r(i)%release()
331 : END DO
332 222 : DEALLOCATE (tau_r)
333 : END IF
334 3196 : ALLOCATE (tau_r(nspins))
335 780 : CALL qs_rho_set(rho, tau_r=tau_r)
336 1636 : DO i = 1, nspins
337 1636 : CALL auxbas_pw_pool%create_pw(tau_r(i))
338 : END DO
339 : END IF
340 :
341 : ! tau_g
342 780 : IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(tau_g)) THEN
343 780 : IF (ASSOCIATED(tau_g)) THEN
344 458 : DO i = 1, SIZE(tau_g)
345 458 : CALL tau_g(i)%release()
346 : END DO
347 222 : DEALLOCATE (tau_g)
348 : END IF
349 3196 : ALLOCATE (tau_g(nspins))
350 780 : CALL qs_rho_set(rho, tau_g=tau_g)
351 1636 : DO i = 1, nspins
352 1636 : CALL auxbas_pw_pool%create_pw(tau_g(i))
353 : END DO
354 : END IF
355 : END IF ! use_kinetic_energy_density
356 :
357 35554 : CALL timestop(handle)
358 :
359 35554 : END SUBROUTINE qs_rho_rebuild
360 :
361 : ! **************************************************************************************************
362 : !> \brief updates rho_r and rho_g to the rho%rho_ao.
363 : !> if use_kinetic_energy_density also computes tau_r and tau_g
364 : !> this works for all ground state and ground state response methods
365 : !> \param rho_struct the rho structure that should be updated
366 : !> \param qs_env the qs_env rho_struct refers to
367 : !> the integrated charge in r space
368 : !> \param rho_xc_external ...
369 : !> \param local_rho_set ...
370 : !> \param task_list_external external task list
371 : !> \param task_list_external_soft external task list (soft_version)
372 : !> \param pw_env_external external plane wave environment
373 : !> \param para_env_external external MPI environment
374 : !> \par History
375 : !> 08.2002 created [fawzi]
376 : !> \author Fawzi Mohamed
377 : ! **************************************************************************************************
378 303067 : SUBROUTINE qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, &
379 : task_list_external, task_list_external_soft, &
380 : pw_env_external, para_env_external)
381 : TYPE(qs_rho_type), INTENT(INOUT) :: rho_struct
382 : TYPE(qs_environment_type), POINTER :: qs_env
383 : TYPE(qs_rho_type), OPTIONAL, POINTER :: rho_xc_external
384 : TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set
385 : TYPE(task_list_type), OPTIONAL, POINTER :: task_list_external, &
386 : task_list_external_soft
387 : TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
388 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_external
389 :
390 303067 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
391 303067 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
392 303067 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
393 303067 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
394 : TYPE(dft_control_type), POINTER :: dft_control
395 : TYPE(harris_type), POINTER :: harris_env
396 : TYPE(kpoint_type), POINTER :: kpoints
397 : TYPE(lri_density_type), POINTER :: lri_density
398 : TYPE(lri_environment_type), POINTER :: lri_env
399 : TYPE(mp_para_env_type), POINTER :: para_env
400 : TYPE(qs_ks_env_type), POINTER :: ks_env
401 :
402 : CALL get_qs_env(qs_env, dft_control=dft_control, &
403 : atomic_kind_set=atomic_kind_set, &
404 303067 : para_env=para_env)
405 303067 : IF (PRESENT(para_env_external)) para_env => para_env_external
406 :
407 303067 : IF (qs_env%harris_method) THEN
408 96 : CALL get_qs_env(qs_env, harris_env=harris_env)
409 96 : CALL calculate_harris_density(qs_env, harris_env, rho_struct)
410 96 : CALL qs_rho_set(rho_struct, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
411 :
412 : ELSE IF (dft_control%qs_control%semi_empirical .OR. &
413 302971 : dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
414 :
415 135858 : CALL qs_rho_set(rho_struct, rho_r_valid=.FALSE., rho_g_valid=.FALSE.)
416 :
417 167113 : ELSE IF (dft_control%qs_control%lrigpw) THEN
418 672 : CPASSERT(.NOT. dft_control%use_kinetic_energy_density)
419 672 : CPASSERT(.NOT. dft_control%drho_by_collocation)
420 672 : CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
421 672 : CALL get_qs_env(qs_env, ks_env=ks_env)
422 672 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
423 672 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
424 672 : CALL get_qs_env(qs_env, lri_env=lri_env, lri_density=lri_density)
425 : CALL calculate_lri_densities(lri_env, lri_density, qs_env, rho_ao_kp, cell_to_index, &
426 : lri_rho_struct=rho_struct, &
427 : atomic_kind_set=atomic_kind_set, &
428 : para_env=para_env, &
429 672 : response_density=.FALSE.)
430 672 : CALL set_qs_env(qs_env, lri_density=lri_density)
431 672 : CALL qs_rho_set(rho_struct, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
432 :
433 166441 : ELSE IF (dft_control%qs_control%rigpw) THEN
434 26 : CPASSERT(.NOT. dft_control%use_kinetic_energy_density)
435 26 : CPASSERT(.NOT. dft_control%drho_by_collocation)
436 26 : CALL get_qs_env(qs_env, lri_env=lri_env)
437 26 : CALL qs_rho_get(rho_struct, rho_ao=rho_ao)
438 : CALL calculate_ri_densities(lri_env, qs_env, rho_ao, &
439 : lri_rho_struct=rho_struct, &
440 : atomic_kind_set=atomic_kind_set, &
441 26 : para_env=para_env)
442 26 : CALL qs_rho_set(rho_struct, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
443 :
444 : ELSE
445 : CALL qs_rho_update_rho_low(rho_struct=rho_struct, qs_env=qs_env, &
446 : rho_xc_external=rho_xc_external, &
447 : local_rho_set=local_rho_set, &
448 : task_list_external=task_list_external, &
449 : task_list_external_soft=task_list_external_soft, &
450 : pw_env_external=pw_env_external, &
451 166415 : para_env_external=para_env_external)
452 :
453 : END IF
454 :
455 303067 : END SUBROUTINE qs_rho_update_rho
456 :
457 : ! **************************************************************************************************
458 : !> \brief updates rho_r and rho_g to the rho%rho_ao.
459 : !> if use_kinetic_energy_density also computes tau_r and tau_g
460 : !> \param rho_struct the rho structure that should be updated
461 : !> \param qs_env the qs_env rho_struct refers to
462 : !> the integrated charge in r space
463 : !> \param rho_xc_external rho structure for GAPW_XC
464 : !> \param local_rho_set ...
465 : !> \param pw_env_external external plane wave environment
466 : !> \param task_list_external external task list (use for default and GAPW)
467 : !> \param task_list_external_soft external task list (soft density for GAPW_XC)
468 : !> \param para_env_external ...
469 : !> \par History
470 : !> 08.2002 created [fawzi]
471 : !> \author Fawzi Mohamed
472 : ! **************************************************************************************************
473 166415 : SUBROUTINE qs_rho_update_rho_low(rho_struct, qs_env, rho_xc_external, &
474 : local_rho_set, pw_env_external, &
475 : task_list_external, task_list_external_soft, &
476 : para_env_external)
477 : TYPE(qs_rho_type), INTENT(INOUT) :: rho_struct
478 : TYPE(qs_environment_type), POINTER :: qs_env
479 : TYPE(qs_rho_type), OPTIONAL, POINTER :: rho_xc_external
480 : TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set
481 : TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
482 : TYPE(task_list_type), OPTIONAL, POINTER :: task_list_external, &
483 : task_list_external_soft
484 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_external
485 :
486 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_rho_update_rho_low'
487 :
488 : INTEGER :: handle, img, ispin, nimg, nspins
489 : LOGICAL :: gapw, gapw_xc
490 : REAL(KIND=dp) :: dum
491 166415 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r, tot_rho_r_xc
492 166415 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
493 166415 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
494 166415 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp, rho_xc_ao
495 : TYPE(dft_control_type), POINTER :: dft_control
496 : TYPE(mp_para_env_type), POINTER :: para_env
497 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
498 166415 : POINTER :: sab
499 : TYPE(oce_matrix_type), POINTER :: oce
500 166415 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rho_xc_g, tau_g, tau_xc_g
501 166415 : TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER :: drho_g, drho_xc_g
502 : TYPE(pw_env_type), POINTER :: pw_env
503 166415 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, rho_xc_r, tau_r, tau_xc_r
504 166415 : TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER :: drho_r, drho_xc_r
505 166415 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
506 : TYPE(qs_ks_env_type), POINTER :: ks_env
507 : TYPE(qs_rho_type), POINTER :: rho_xc
508 166415 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
509 : TYPE(task_list_type), POINTER :: task_list
510 :
511 166415 : CALL timeset(routineN, handle)
512 :
513 166415 : NULLIFY (dft_control, rho_xc, ks_env, rho_ao, rho_r, rho_g, drho_r, drho_g, tau_r, tau_g)
514 166415 : 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)
515 166415 : NULLIFY (para_env, pw_env, atomic_kind_set)
516 :
517 : CALL get_qs_env(qs_env, &
518 : ks_env=ks_env, &
519 : dft_control=dft_control, &
520 166415 : atomic_kind_set=atomic_kind_set)
521 :
522 : CALL qs_rho_get(rho_struct, &
523 : rho_r=rho_r, &
524 : rho_g=rho_g, &
525 : tot_rho_r=tot_rho_r, &
526 : drho_r=drho_r, &
527 : drho_g=drho_g, &
528 : tau_r=tau_r, &
529 166415 : tau_g=tau_g)
530 :
531 : CALL get_qs_env(qs_env, task_list=task_list, &
532 166415 : para_env=para_env, pw_env=pw_env)
533 166415 : IF (PRESENT(pw_env_external)) pw_env => pw_env_external
534 166415 : IF (PRESENT(task_list_external)) task_list => task_list_external
535 166415 : IF (PRESENT(para_env_external)) para_env => para_env_external
536 :
537 166415 : nspins = dft_control%nspins
538 166415 : nimg = dft_control%nimages
539 166415 : gapw = dft_control%qs_control%gapw
540 166415 : gapw_xc = dft_control%qs_control%gapw_xc
541 :
542 166415 : CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
543 365413 : DO ispin = 1, nspins
544 198998 : rho_ao => rho_ao_kp(ispin, :)
545 : CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
546 : rho=rho_r(ispin), &
547 : rho_gspace=rho_g(ispin), &
548 : total_rho=tot_rho_r(ispin), &
549 : ks_env=ks_env, soft_valid=gapw, &
550 : task_list_external=task_list_external, &
551 365413 : pw_env_external=pw_env_external)
552 : END DO
553 166415 : CALL qs_rho_set(rho_struct, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
554 :
555 166415 : IF (gapw_xc) THEN
556 6180 : IF (PRESENT(rho_xc_external)) THEN
557 1138 : rho_xc => rho_xc_external
558 : ELSE
559 5042 : CALL get_qs_env(qs_env=qs_env, rho_xc=rho_xc)
560 : END IF
561 : CALL qs_rho_get(rho_xc, &
562 : rho_ao_kp=rho_xc_ao, &
563 : rho_r=rho_xc_r, &
564 : rho_g=rho_xc_g, &
565 6180 : tot_rho_r=tot_rho_r_xc)
566 : ! copy rho_ao into rho_xc_ao
567 12818 : DO ispin = 1, nspins
568 29258 : DO img = 1, nimg
569 23078 : CALL dbcsr_copy(rho_xc_ao(ispin, img)%matrix, rho_ao_kp(ispin, img)%matrix)
570 : END DO
571 : END DO
572 12818 : DO ispin = 1, nspins
573 6638 : rho_ao => rho_xc_ao(ispin, :)
574 : CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
575 : rho=rho_xc_r(ispin), &
576 : rho_gspace=rho_xc_g(ispin), &
577 : total_rho=tot_rho_r_xc(ispin), &
578 : ks_env=ks_env, soft_valid=gapw_xc, &
579 : task_list_external=task_list_external_soft, &
580 12818 : pw_env_external=pw_env_external)
581 : END DO
582 6180 : CALL qs_rho_set(rho_xc, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
583 : END IF
584 :
585 : ! GAPW o GAPW_XC require the calculation of hard and soft local densities
586 166415 : IF (gapw .OR. gapw_xc) THEN
587 : CALL get_qs_env(qs_env=qs_env, &
588 : rho_atom_set=rho_atom_set, &
589 : qs_kind_set=qs_kind_set, &
590 36626 : oce=oce, sab_orb=sab)
591 36626 : IF (PRESENT(local_rho_set)) rho_atom_set => local_rho_set%rho_atom_set
592 36626 : CPASSERT(ASSOCIATED(rho_atom_set))
593 36626 : CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
594 36626 : CALL calculate_rho_atom_coeff(qs_env, rho_ao_kp, rho_atom_set, qs_kind_set, oce, sab, para_env)
595 : END IF
596 :
597 166415 : IF (.NOT. gapw_xc) THEN
598 : ! if needed compute also the gradient of the density
599 160235 : IF (dft_control%drho_by_collocation) THEN
600 0 : CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
601 0 : CPASSERT(.NOT. PRESENT(task_list_external))
602 0 : DO ispin = 1, nspins
603 0 : rho_ao => rho_ao_kp(ispin, :)
604 : CALL calculate_drho_elec(matrix_p_kp=rho_ao, &
605 : drho=drho_r(:, ispin), &
606 : drho_gspace=drho_g(:, ispin), &
607 0 : qs_env=qs_env, soft_valid=gapw)
608 : END DO
609 0 : CALL qs_rho_set(rho_struct, drho_r_valid=.TRUE., drho_g_valid=.TRUE.)
610 : END IF
611 : ! if needed compute also the kinetic energy density
612 160235 : IF (dft_control%use_kinetic_energy_density) THEN
613 5962 : CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
614 12702 : DO ispin = 1, nspins
615 6740 : rho_ao => rho_ao_kp(ispin, :)
616 : CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
617 : rho=tau_r(ispin), &
618 : rho_gspace=tau_g(ispin), &
619 : total_rho=dum, & ! presumably not meaningful
620 : ks_env=ks_env, soft_valid=gapw, &
621 : compute_tau=.TRUE., &
622 : task_list_external=task_list_external, &
623 12702 : pw_env_external=pw_env_external)
624 : END DO
625 5962 : CALL qs_rho_set(rho_struct, tau_r_valid=.TRUE., tau_g_valid=.TRUE.)
626 : END IF
627 : ELSE
628 : CALL qs_rho_get(rho_xc, &
629 : drho_r=drho_xc_r, &
630 : drho_g=drho_xc_g, &
631 : tau_r=tau_xc_r, &
632 6180 : tau_g=tau_xc_g)
633 : ! if needed compute also the gradient of the density
634 6180 : IF (dft_control%drho_by_collocation) THEN
635 0 : CPASSERT(.NOT. PRESENT(task_list_external))
636 0 : DO ispin = 1, nspins
637 0 : rho_ao => rho_xc_ao(ispin, :)
638 : CALL calculate_drho_elec(matrix_p_kp=rho_ao, &
639 : drho=drho_xc_r(:, ispin), &
640 : drho_gspace=drho_xc_g(:, ispin), &
641 0 : qs_env=qs_env, soft_valid=gapw_xc)
642 : END DO
643 0 : CALL qs_rho_set(rho_xc, drho_r_valid=.TRUE., drho_g_valid=.TRUE.)
644 : END IF
645 : ! if needed compute also the kinetic energy density
646 6180 : IF (dft_control%use_kinetic_energy_density) THEN
647 724 : DO ispin = 1, nspins
648 362 : rho_ao => rho_xc_ao(ispin, :)
649 : CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
650 : rho=tau_xc_r(ispin), &
651 : rho_gspace=tau_xc_g(ispin), &
652 : ks_env=ks_env, soft_valid=gapw_xc, &
653 : compute_tau=.TRUE., &
654 : task_list_external=task_list_external_soft, &
655 724 : pw_env_external=pw_env_external)
656 : END DO
657 362 : CALL qs_rho_set(rho_xc, tau_r_valid=.TRUE., tau_g_valid=.TRUE.)
658 : END IF
659 : END IF
660 :
661 166415 : CALL timestop(handle)
662 :
663 166415 : END SUBROUTINE qs_rho_update_rho_low
664 :
665 : ! **************************************************************************************************
666 : !> \brief updates rho_r and rho_g to the rho%rho_ao.
667 : !> if use_kinetic_energy_density also computes tau_r and tau_g
668 : !> \param rho_struct the rho structure that should be updated
669 : !> \param qs_env the qs_env rho_struct refers to
670 : !> the integrated charge in r space
671 : !> \param pw_env_external external plane wave environment
672 : !> \param task_list_external external task list
673 : !> \param para_env_external ...
674 : !> \param tddfpt_lri_env ...
675 : !> \param tddfpt_lri_density ...
676 : !> \par History
677 : !> 08.2002 created [fawzi]
678 : !> \author Fawzi Mohamed
679 : ! **************************************************************************************************
680 172 : SUBROUTINE qs_rho_update_tddfpt(rho_struct, qs_env, pw_env_external, task_list_external, &
681 : para_env_external, tddfpt_lri_env, tddfpt_lri_density)
682 : TYPE(qs_rho_type), INTENT(INOUT) :: rho_struct
683 : TYPE(qs_environment_type), POINTER :: qs_env
684 : TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
685 : TYPE(task_list_type), OPTIONAL, POINTER :: task_list_external
686 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env_external
687 : TYPE(lri_environment_type), OPTIONAL, POINTER :: tddfpt_lri_env
688 : TYPE(lri_density_type), OPTIONAL, POINTER :: tddfpt_lri_density
689 :
690 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_rho_update_tddfpt'
691 :
692 : INTEGER :: handle, ispin, nspins
693 172 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
694 : LOGICAL :: lri_response
695 172 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
696 172 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
697 172 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
698 172 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
699 : TYPE(dft_control_type), POINTER :: dft_control
700 : TYPE(kpoint_type), POINTER :: kpoints
701 : TYPE(mp_para_env_type), POINTER :: para_env
702 172 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
703 : TYPE(pw_env_type), POINTER :: pw_env
704 172 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
705 : TYPE(qs_ks_env_type), POINTER :: ks_env
706 : TYPE(task_list_type), POINTER :: task_list
707 :
708 172 : CALL timeset(routineN, handle)
709 :
710 : CALL get_qs_env(qs_env, &
711 : ks_env=ks_env, &
712 : dft_control=dft_control, &
713 : atomic_kind_set=atomic_kind_set, &
714 : task_list=task_list, &
715 : para_env=para_env, &
716 172 : pw_env=pw_env)
717 172 : IF (PRESENT(pw_env_external)) pw_env => pw_env_external
718 172 : IF (PRESENT(task_list_external)) task_list => task_list_external
719 172 : IF (PRESENT(para_env_external)) para_env => para_env_external
720 :
721 : CALL qs_rho_get(rho_struct, &
722 : rho_r=rho_r, &
723 : rho_g=rho_g, &
724 172 : tot_rho_r=tot_rho_r)
725 :
726 172 : nspins = dft_control%nspins
727 :
728 172 : lri_response = PRESENT(tddfpt_lri_env)
729 172 : IF (lri_response) THEN
730 172 : CPASSERT(PRESENT(tddfpt_lri_density))
731 : END IF
732 :
733 172 : CPASSERT(.NOT. dft_control%drho_by_collocation)
734 172 : CPASSERT(.NOT. dft_control%use_kinetic_energy_density)
735 172 : CPASSERT(.NOT. dft_control%qs_control%gapw)
736 172 : CPASSERT(.NOT. dft_control%qs_control%gapw_xc)
737 :
738 172 : IF (lri_response) THEN
739 172 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
740 172 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
741 172 : CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
742 : CALL calculate_lri_densities(tddfpt_lri_env, tddfpt_lri_density, qs_env, rho_ao_kp, cell_to_index, &
743 : lri_rho_struct=rho_struct, &
744 : atomic_kind_set=atomic_kind_set, &
745 : para_env=para_env, &
746 172 : response_density=lri_response)
747 172 : CALL qs_rho_set(rho_struct, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
748 : ELSE
749 0 : CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
750 0 : DO ispin = 1, nspins
751 0 : rho_ao => rho_ao_kp(ispin, :)
752 : CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
753 : rho=rho_r(ispin), &
754 : rho_gspace=rho_g(ispin), &
755 : total_rho=tot_rho_r(ispin), &
756 : ks_env=ks_env, &
757 : task_list_external=task_list_external, &
758 0 : pw_env_external=pw_env_external)
759 : END DO
760 0 : CALL qs_rho_set(rho_struct, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
761 : END IF
762 :
763 172 : CALL timestop(handle)
764 :
765 172 : END SUBROUTINE qs_rho_update_tddfpt
766 :
767 : ! **************************************************************************************************
768 : !> \brief Allocate a density structure and fill it with data from an input structure
769 : !> SIZE(rho_input) == mspin == 1 direct copy
770 : !> SIZE(rho_input) == mspin == 2 direct copy of alpha and beta spin
771 : !> SIZE(rho_input) == 1 AND mspin == 2 copy rho/2 into alpha and beta spin
772 : !> \param rho_input ...
773 : !> \param rho_output ...
774 : !> \param auxbas_pw_pool ...
775 : !> \param mspin ...
776 : !> \param factor ...
777 : ! **************************************************************************************************
778 35252 : SUBROUTINE qs_rho_copy(rho_input, rho_output, auxbas_pw_pool, mspin, factor)
779 :
780 : TYPE(qs_rho_type), INTENT(IN) :: rho_input
781 : TYPE(qs_rho_type), INTENT(INOUT) :: rho_output
782 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
783 : INTEGER, INTENT(IN) :: mspin
784 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: factor
785 :
786 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_rho_copy'
787 :
788 : INTEGER :: handle, i, j, nspins
789 : LOGICAL :: complex_rho_ao, drho_g_valid_in, drho_r_valid_in, rho_g_valid_in, rho_r_valid_in, &
790 : soft_valid_in, tau_g_valid_in, tau_r_valid_in
791 : REAL(KIND=dp) :: ospin
792 17626 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_g_in, tot_rho_g_out, &
793 17626 : tot_rho_r_in, tot_rho_r_out
794 17626 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_im_in, rho_ao_im_out, rho_ao_in, &
795 17626 : rho_ao_out
796 17626 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp_in
797 17626 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_in, rho_g_out, tau_g_in, tau_g_out
798 17626 : TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER :: drho_g_in, drho_g_out
799 17626 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_in, rho_r_out, tau_r_in, tau_r_out
800 17626 : TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER :: drho_r_in, drho_r_out
801 : TYPE(pw_r3d_rs_type), POINTER :: rho_r_sccs_in, rho_r_sccs_out
802 :
803 17626 : CALL timeset(routineN, handle)
804 :
805 17626 : CPASSERT(mspin == 1 .OR. mspin == 2)
806 17626 : ospin = 1._dp/REAL(mspin, KIND=dp)
807 17626 : IF (PRESENT(factor)) THEN
808 774 : ospin = ospin*factor
809 : END IF
810 :
811 17626 : CALL qs_rho_clear(rho_output)
812 :
813 17626 : NULLIFY (rho_ao_in, rho_ao_kp_in, rho_ao_im_in, rho_r_in, rho_g_in, drho_r_in, &
814 17626 : drho_g_in, tau_r_in, tau_g_in, tot_rho_r_in, tot_rho_g_in, rho_r_sccs_in)
815 :
816 : CALL qs_rho_get(rho_input, &
817 : rho_ao=rho_ao_in, &
818 : rho_ao_kp=rho_ao_kp_in, &
819 : rho_ao_im=rho_ao_im_in, &
820 : rho_r=rho_r_in, &
821 : rho_g=rho_g_in, &
822 : drho_r=drho_r_in, &
823 : drho_g=drho_g_in, &
824 : tau_r=tau_r_in, &
825 : tau_g=tau_g_in, &
826 : tot_rho_r=tot_rho_r_in, &
827 : tot_rho_g=tot_rho_g_in, &
828 : rho_g_valid=rho_g_valid_in, &
829 : rho_r_valid=rho_r_valid_in, &
830 : drho_g_valid=drho_g_valid_in, &
831 : drho_r_valid=drho_r_valid_in, &
832 : tau_r_valid=tau_r_valid_in, &
833 : tau_g_valid=tau_g_valid_in, &
834 : rho_r_sccs=rho_r_sccs_in, &
835 : soft_valid=soft_valid_in, &
836 17626 : complex_rho_ao=complex_rho_ao)
837 :
838 17626 : NULLIFY (rho_ao_out, rho_ao_im_out, rho_r_out, rho_g_out, drho_r_out, &
839 17626 : drho_g_out, tau_r_out, tau_g_out, tot_rho_r_out, tot_rho_g_out, rho_r_sccs_out)
840 : ! rho_ao
841 17626 : IF (ASSOCIATED(rho_ao_in)) THEN
842 17626 : nspins = SIZE(rho_ao_in)
843 17626 : CPASSERT(mspin >= nspins)
844 17626 : CALL dbcsr_allocate_matrix_set(rho_ao_out, mspin)
845 17626 : CALL qs_rho_set(rho_output, rho_ao=rho_ao_out)
846 17626 : IF (mspin > nspins) THEN
847 5850 : DO i = 1, mspin
848 3900 : ALLOCATE (rho_ao_out(i)%matrix)
849 3900 : CALL dbcsr_copy(rho_ao_out(i)%matrix, rho_ao_in(1)%matrix, name="RHO copy")
850 5850 : CALL dbcsr_scale(rho_ao_out(i)%matrix, ospin)
851 : END DO
852 : ELSE
853 33312 : DO i = 1, nspins
854 17636 : ALLOCATE (rho_ao_out(i)%matrix)
855 33312 : CALL dbcsr_copy(rho_ao_out(i)%matrix, rho_ao_in(i)%matrix, name="RHO copy")
856 : END DO
857 : END IF
858 : END IF
859 :
860 : ! rho_ao_kp
861 : ! only for non-kp, we could probably just copy this pointer, should work also for non-kp?
862 : !IF (ASSOCIATED(rho_ao_kp_in)) THEN
863 : ! CPABORT("Copy not available")
864 : !END IF
865 :
866 : ! rho_ao_im
867 17626 : IF (ASSOCIATED(rho_ao_im_in)) THEN
868 0 : nspins = SIZE(rho_ao_im_in)
869 0 : CPASSERT(mspin >= nspins)
870 0 : CALL dbcsr_allocate_matrix_set(rho_ao_im_out, mspin)
871 0 : CALL qs_rho_set(rho_output, rho_ao_im=rho_ao_im_out)
872 0 : IF (mspin > nspins) THEN
873 0 : DO i = 1, mspin
874 0 : ALLOCATE (rho_ao_im_out(i)%matrix)
875 0 : CALL dbcsr_copy(rho_ao_im_out(i)%matrix, rho_ao_im_in(1)%matrix, name="RHO copy")
876 0 : CALL dbcsr_scale(rho_ao_im_out(i)%matrix, ospin)
877 : END DO
878 : ELSE
879 0 : DO i = 1, nspins
880 0 : ALLOCATE (rho_ao_im_out(i)%matrix)
881 0 : CALL dbcsr_copy(rho_ao_im_out(i)%matrix, rho_ao_im_in(i)%matrix, name="RHO copy")
882 : END DO
883 : END IF
884 : END IF
885 :
886 : ! rho_r
887 17626 : IF (ASSOCIATED(rho_r_in)) THEN
888 17626 : nspins = SIZE(rho_r_in)
889 17626 : CPASSERT(mspin >= nspins)
890 74414 : ALLOCATE (rho_r_out(mspin))
891 17626 : CALL qs_rho_set(rho_output, rho_r=rho_r_out)
892 17626 : IF (mspin > nspins) THEN
893 5850 : DO i = 1, mspin
894 3900 : CALL auxbas_pw_pool%create_pw(rho_r_out(i))
895 3900 : CALL pw_copy(rho_r_in(1), rho_r_out(i))
896 5850 : CALL pw_scale(rho_r_out(i), ospin)
897 : END DO
898 : ELSE
899 33312 : DO i = 1, nspins
900 17636 : CALL auxbas_pw_pool%create_pw(rho_r_out(i))
901 33312 : CALL pw_copy(rho_r_in(i), rho_r_out(i))
902 : END DO
903 : END IF
904 : END IF
905 :
906 : ! rho_g
907 17626 : IF (ASSOCIATED(rho_g_in)) THEN
908 17626 : nspins = SIZE(rho_g_in)
909 17626 : CPASSERT(mspin >= nspins)
910 74414 : ALLOCATE (rho_g_out(mspin))
911 17626 : CALL qs_rho_set(rho_output, rho_g=rho_g_out)
912 17626 : IF (mspin > nspins) THEN
913 5850 : DO i = 1, mspin
914 3900 : CALL auxbas_pw_pool%create_pw(rho_g_out(i))
915 3900 : CALL pw_copy(rho_g_in(1), rho_g_out(i))
916 5850 : CALL pw_scale(rho_g_out(i), ospin)
917 : END DO
918 : ELSE
919 33312 : DO i = 1, nspins
920 17636 : CALL auxbas_pw_pool%create_pw(rho_g_out(i))
921 33312 : CALL pw_copy(rho_g_in(i), rho_g_out(i))
922 : END DO
923 : END IF
924 : END IF
925 :
926 : ! SCCS
927 17626 : IF (ASSOCIATED(rho_r_sccs_in)) THEN
928 0 : CALL qs_rho_set(rho_output, rho_r_sccs=rho_r_sccs_out)
929 0 : CALL auxbas_pw_pool%create_pw(rho_r_sccs_out)
930 0 : CALL pw_copy(rho_r_sccs_in, rho_r_sccs_out)
931 : END IF
932 :
933 : ! drho_r
934 17626 : IF (ASSOCIATED(drho_r_in)) THEN
935 0 : nspins = SIZE(drho_r_in)
936 0 : CPASSERT(mspin >= nspins)
937 0 : ALLOCATE (drho_r_out(3, mspin))
938 0 : CALL qs_rho_set(rho_output, drho_r=drho_r_out)
939 0 : IF (mspin > nspins) THEN
940 0 : DO j = 1, mspin
941 0 : DO i = 1, 3
942 0 : CALL auxbas_pw_pool%create_pw(drho_r_out(i, j))
943 0 : CALL pw_copy(drho_r_in(i, 1), drho_r_out(i, j))
944 0 : CALL pw_scale(drho_r_out(i, j), ospin)
945 : END DO
946 : END DO
947 : ELSE
948 0 : DO j = 1, nspins
949 0 : DO i = 1, 3
950 0 : CALL auxbas_pw_pool%create_pw(drho_r_out(i, j))
951 0 : CALL pw_copy(drho_r_in(i, j), drho_r_out(i, j))
952 : END DO
953 : END DO
954 : END IF
955 : END IF
956 :
957 : ! drho_g
958 17626 : IF (ASSOCIATED(drho_g_in)) THEN
959 0 : nspins = SIZE(drho_g_in)
960 0 : CPASSERT(mspin >= nspins)
961 0 : ALLOCATE (drho_g_out(3, mspin))
962 0 : CALL qs_rho_set(rho_output, drho_g=drho_g_out)
963 0 : IF (mspin > nspins) THEN
964 0 : DO j = 1, mspin
965 0 : DO i = 1, 3
966 0 : CALL auxbas_pw_pool%create_pw(drho_g_out(i, j))
967 0 : CALL pw_copy(drho_g_in(i, 1), drho_g_out(i, j))
968 0 : CALL pw_scale(drho_g_out(i, j), ospin)
969 : END DO
970 : END DO
971 : ELSE
972 0 : DO j = 1, nspins
973 0 : DO i = 1, 3
974 0 : CALL auxbas_pw_pool%create_pw(drho_g_out(i, j))
975 0 : CALL pw_copy(drho_g_in(i, j), drho_g_out(i, j))
976 : END DO
977 : END DO
978 : END IF
979 : END IF
980 :
981 : ! tau_r
982 17626 : IF (ASSOCIATED(tau_r_in)) THEN
983 6 : nspins = SIZE(tau_r_in)
984 6 : CPASSERT(mspin >= nspins)
985 24 : ALLOCATE (tau_r_out(mspin))
986 6 : CALL qs_rho_set(rho_output, tau_r=tau_r_out)
987 6 : IF (mspin > nspins) THEN
988 0 : DO i = 1, mspin
989 0 : CALL auxbas_pw_pool%create_pw(tau_r_out(i))
990 0 : CALL pw_copy(tau_r_in(1), tau_r_out(i))
991 0 : CALL pw_scale(tau_r_out(i), ospin)
992 : END DO
993 : ELSE
994 12 : DO i = 1, nspins
995 6 : CALL auxbas_pw_pool%create_pw(tau_r_out(i))
996 12 : CALL pw_copy(tau_r_in(i), tau_r_out(i))
997 : END DO
998 : END IF
999 : END IF
1000 :
1001 : ! tau_g
1002 17626 : IF (ASSOCIATED(tau_g_in)) THEN
1003 6 : nspins = SIZE(tau_g_in)
1004 6 : CPASSERT(mspin >= nspins)
1005 24 : ALLOCATE (tau_g_out(mspin))
1006 6 : CALL qs_rho_set(rho_output, tau_g=tau_g_out)
1007 6 : IF (mspin > nspins) THEN
1008 0 : DO i = 1, mspin
1009 0 : CALL auxbas_pw_pool%create_pw(tau_g_out(i))
1010 0 : CALL pw_copy(tau_g_in(1), tau_g_out(i))
1011 0 : CALL pw_scale(tau_g_out(i), ospin)
1012 : END DO
1013 : ELSE
1014 12 : DO i = 1, nspins
1015 6 : CALL auxbas_pw_pool%create_pw(tau_g_out(i))
1016 12 : CALL pw_copy(tau_g_in(i), tau_g_out(i))
1017 : END DO
1018 : END IF
1019 : END IF
1020 :
1021 : ! tot_rho_r
1022 17626 : IF (ASSOCIATED(tot_rho_r_in)) THEN
1023 17620 : nspins = SIZE(tot_rho_r_in)
1024 17620 : CPASSERT(mspin >= nspins)
1025 52860 : ALLOCATE (tot_rho_r_out(mspin))
1026 17620 : CALL qs_rho_set(rho_output, tot_rho_r=tot_rho_r_out)
1027 17620 : IF (mspin > nspins) THEN
1028 5832 : DO i = 1, mspin
1029 5832 : tot_rho_r_out(i) = tot_rho_r_in(1)*ospin
1030 : END DO
1031 : ELSE
1032 33312 : DO i = 1, nspins
1033 33312 : tot_rho_r_out(i) = tot_rho_r_in(i)
1034 : END DO
1035 : END IF
1036 : END IF
1037 :
1038 : ! tot_rho_g
1039 17626 : IF (ASSOCIATED(tot_rho_g_in)) THEN
1040 0 : nspins = SIZE(tot_rho_g_in)
1041 0 : CPASSERT(mspin >= nspins)
1042 0 : ALLOCATE (tot_rho_g_out(mspin))
1043 0 : CALL qs_rho_set(rho_output, tot_rho_g=tot_rho_g_out)
1044 0 : IF (mspin > nspins) THEN
1045 0 : DO i = 1, mspin
1046 0 : tot_rho_g_out(i) = tot_rho_g_in(1)*ospin
1047 : END DO
1048 : ELSE
1049 0 : DO i = 1, nspins
1050 0 : tot_rho_g_out(i) = tot_rho_g_in(i)
1051 : END DO
1052 : END IF
1053 : END IF
1054 :
1055 : CALL qs_rho_set(rho_output, &
1056 : rho_g_valid=rho_g_valid_in, &
1057 : rho_r_valid=rho_r_valid_in, &
1058 : drho_g_valid=drho_g_valid_in, &
1059 : drho_r_valid=drho_r_valid_in, &
1060 : tau_r_valid=tau_r_valid_in, &
1061 : tau_g_valid=tau_g_valid_in, &
1062 : soft_valid=soft_valid_in, &
1063 17626 : complex_rho_ao=complex_rho_ao)
1064 :
1065 17626 : CALL timestop(handle)
1066 :
1067 17626 : END SUBROUTINE qs_rho_copy
1068 :
1069 : ! **************************************************************************************************
1070 : !> \brief Allocate a density structure and fill it with data from an input structure
1071 : !> Transfer all data to input pw_pool
1072 : !> \param rho_input ...
1073 : !> \param rho_output ...
1074 : !> \param in_pw_pool ...
1075 : !> \param out_pw_pool ...
1076 : ! **************************************************************************************************
1077 552 : SUBROUTINE qs_rho_transfer(rho_input, rho_output, in_pw_pool, out_pw_pool)
1078 :
1079 : TYPE(qs_rho_type), INTENT(IN) :: rho_input
1080 : TYPE(qs_rho_type), INTENT(INOUT) :: rho_output
1081 : TYPE(pw_pool_type), POINTER :: in_pw_pool, out_pw_pool
1082 :
1083 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_rho_transfer'
1084 :
1085 : INTEGER :: handle, i, j, nspins
1086 : LOGICAL :: complex_rho_ao, drho_g_valid_in, drho_r_valid_in, rho_g_valid_in, rho_r_valid_in, &
1087 : soft_valid_in, tau_g_valid_in, tau_r_valid_in
1088 276 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_g_in, tot_rho_g_out, &
1089 276 : tot_rho_r_in, tot_rho_r_out
1090 276 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_im_in, rho_ao_im_out, rho_ao_in, &
1091 276 : rho_ao_out
1092 276 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp_in
1093 276 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_in, rho_g_out, tau_g_in, tau_g_out
1094 276 : TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER :: drho_g_in, drho_g_out
1095 276 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_in, rho_r_out, tau_r_in, tau_r_out
1096 276 : TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER :: drho_r_in, drho_r_out
1097 : TYPE(pw_r3d_rs_type), POINTER :: rho_r_sccs_in, rho_r_sccs_out
1098 :
1099 276 : CALL timeset(routineN, handle)
1100 :
1101 276 : CALL qs_rho_clear(rho_output)
1102 :
1103 276 : NULLIFY (rho_ao_in, rho_ao_kp_in, rho_ao_im_in, rho_r_in, rho_g_in, drho_r_in, &
1104 276 : drho_g_in, tau_r_in, tau_g_in, tot_rho_r_in, tot_rho_g_in, rho_r_sccs_in)
1105 :
1106 : CALL qs_rho_get(rho_input, &
1107 : rho_ao=rho_ao_in, &
1108 : rho_ao_kp=rho_ao_kp_in, &
1109 : rho_ao_im=rho_ao_im_in, &
1110 : rho_r=rho_r_in, &
1111 : rho_g=rho_g_in, &
1112 : drho_r=drho_r_in, &
1113 : drho_g=drho_g_in, &
1114 : tau_r=tau_r_in, &
1115 : tau_g=tau_g_in, &
1116 : tot_rho_r=tot_rho_r_in, &
1117 : tot_rho_g=tot_rho_g_in, &
1118 : rho_g_valid=rho_g_valid_in, &
1119 : rho_r_valid=rho_r_valid_in, &
1120 : drho_g_valid=drho_g_valid_in, &
1121 : drho_r_valid=drho_r_valid_in, &
1122 : tau_r_valid=tau_r_valid_in, &
1123 : tau_g_valid=tau_g_valid_in, &
1124 : rho_r_sccs=rho_r_sccs_in, &
1125 : soft_valid=soft_valid_in, &
1126 276 : complex_rho_ao=complex_rho_ao)
1127 :
1128 276 : NULLIFY (rho_ao_out, rho_ao_im_out, rho_r_out, rho_g_out, drho_r_out, &
1129 276 : drho_g_out, tau_r_out, tau_g_out, tot_rho_r_out, tot_rho_g_out, rho_r_sccs_out)
1130 : ! rho_ao
1131 276 : IF (ASSOCIATED(rho_ao_in)) THEN
1132 260 : nspins = SIZE(rho_ao_in)
1133 260 : CALL dbcsr_allocate_matrix_set(rho_ao_out, nspins)
1134 260 : CALL qs_rho_set(rho_output, rho_ao=rho_ao_out)
1135 520 : DO i = 1, nspins
1136 260 : ALLOCATE (rho_ao_out(i)%matrix)
1137 520 : CALL dbcsr_copy(rho_ao_out(i)%matrix, rho_ao_in(i)%matrix, name="RHO copy")
1138 : END DO
1139 : END IF
1140 :
1141 : ! rho_ao_kp
1142 : ! only for non-kp, we could probably just copy this pointer, should work also for non-kp?
1143 : !IF (ASSOCIATED(rho_ao_kp_in)) THEN
1144 : ! CPABORT("Copy not available")
1145 : !END IF
1146 :
1147 : ! rho_ao_im
1148 276 : IF (ASSOCIATED(rho_ao_im_in)) THEN
1149 0 : nspins = SIZE(rho_ao_im_in)
1150 0 : CALL dbcsr_allocate_matrix_set(rho_ao_im_out, nspins)
1151 0 : CALL qs_rho_set(rho_output, rho_ao_im=rho_ao_im_out)
1152 0 : DO i = 1, nspins
1153 0 : ALLOCATE (rho_ao_im_out(i)%matrix)
1154 0 : CALL dbcsr_copy(rho_ao_im_out(i)%matrix, rho_ao_im_in(i)%matrix, name="RHO copy")
1155 : END DO
1156 : END IF
1157 :
1158 : ! rho_r and rho_g
1159 276 : IF (ASSOCIATED(rho_g_in)) THEN
1160 268 : nspins = SIZE(rho_g_in)
1161 1072 : ALLOCATE (rho_g_out(nspins))
1162 268 : CALL qs_rho_set(rho_output, rho_g=rho_g_out)
1163 268 : IF (ASSOCIATED(rho_r_in)) THEN
1164 1072 : ALLOCATE (rho_r_out(nspins))
1165 268 : CALL qs_rho_set(rho_output, rho_r=rho_r_out)
1166 : END IF
1167 536 : DO i = 1, nspins
1168 268 : CALL out_pw_pool%create_pw(rho_g_out(i))
1169 268 : CALL pw_transfer(rho_g_in(i), rho_g_out(i))
1170 536 : IF (ASSOCIATED(rho_r_in)) THEN
1171 268 : CALL out_pw_pool%create_pw(rho_r_out(i))
1172 268 : CALL pw_transfer(rho_g_out(i), rho_r_out(i))
1173 : END IF
1174 : END DO
1175 8 : ELSE IF (ASSOCIATED(rho_r_in)) THEN
1176 8 : nspins = SIZE(rho_r_in)
1177 32 : ALLOCATE (rho_r_out(nspins))
1178 8 : CALL qs_rho_set(rho_output, rho_r=rho_r_out)
1179 16 : DO i = 1, nspins
1180 8 : CALL out_pw_pool%create_pw(rho_r_out(i))
1181 8 : BLOCK
1182 : TYPE(pw_c1d_gs_type) :: grho_in, grho_out
1183 8 : CALL in_pw_pool%create_pw(grho_in)
1184 8 : CALL out_pw_pool%create_pw(grho_out)
1185 8 : CALL pw_transfer(rho_r_in(i), grho_in)
1186 8 : CALL pw_transfer(grho_in, grho_out)
1187 8 : CALL pw_transfer(grho_out, rho_r_out(i))
1188 8 : CALL out_pw_pool%give_back_pw(grho_out)
1189 16 : CALL in_pw_pool%give_back_pw(grho_in)
1190 : END BLOCK
1191 : END DO
1192 : END IF
1193 :
1194 : ! SCCS
1195 276 : IF (ASSOCIATED(rho_r_sccs_in)) THEN
1196 0 : CALL qs_rho_set(rho_output, rho_r_sccs=rho_r_sccs_out)
1197 0 : CALL out_pw_pool%create_pw(rho_r_sccs_out)
1198 : BLOCK
1199 : TYPE(pw_c1d_gs_type) :: grho_in, grho_out
1200 0 : CALL in_pw_pool%create_pw(grho_in)
1201 0 : CALL out_pw_pool%create_pw(grho_out)
1202 0 : CALL pw_transfer(rho_r_sccs_in, grho_in)
1203 0 : CALL pw_transfer(grho_in, grho_out)
1204 0 : CALL pw_transfer(grho_out, rho_r_sccs_out)
1205 0 : CALL out_pw_pool%give_back_pw(grho_out)
1206 0 : CALL in_pw_pool%give_back_pw(grho_in)
1207 : END BLOCK
1208 : END IF
1209 :
1210 : ! drho_r and drho_g
1211 276 : IF (ASSOCIATED(drho_g_in)) THEN
1212 0 : nspins = SIZE(drho_g_in)
1213 0 : ALLOCATE (drho_g_out(3, nspins))
1214 0 : CALL qs_rho_set(rho_output, drho_g=drho_g_out)
1215 0 : IF (ASSOCIATED(drho_r_in)) THEN
1216 0 : ALLOCATE (drho_r_out(3, nspins))
1217 0 : CALL qs_rho_set(rho_output, drho_r=drho_r_out)
1218 : END IF
1219 0 : DO i = 1, nspins
1220 0 : DO j = 1, 3
1221 0 : CALL out_pw_pool%create_pw(drho_g_out(j, i))
1222 0 : CALL pw_transfer(drho_g_in(j, i), drho_g_out(j, i))
1223 0 : IF (ASSOCIATED(drho_r_in)) THEN
1224 0 : CALL out_pw_pool%create_pw(drho_r_out(j, i))
1225 0 : CALL pw_transfer(drho_g_out(j, i), drho_r_out(j, i))
1226 : END IF
1227 : END DO
1228 : END DO
1229 276 : ELSE IF (ASSOCIATED(drho_r_in)) THEN
1230 0 : nspins = SIZE(drho_r_in)
1231 0 : ALLOCATE (drho_r_out(3, nspins))
1232 0 : CALL qs_rho_set(rho_output, drho_r=drho_r_out)
1233 0 : DO i = 1, nspins
1234 0 : BLOCK
1235 : TYPE(pw_c1d_gs_type) :: grho_in, grho_out
1236 0 : CALL in_pw_pool%create_pw(grho_in)
1237 0 : CALL out_pw_pool%create_pw(grho_out)
1238 0 : DO j = 1, 3
1239 0 : CALL out_pw_pool%create_pw(drho_r_out(j, i))
1240 0 : CALL pw_transfer(drho_r_in(j, i), grho_in)
1241 0 : CALL pw_transfer(grho_in, grho_out)
1242 0 : CALL pw_transfer(grho_out, drho_r_out(j, i))
1243 : END DO
1244 0 : CALL out_pw_pool%give_back_pw(grho_out)
1245 0 : CALL in_pw_pool%give_back_pw(grho_in)
1246 : END BLOCK
1247 : END DO
1248 : END IF
1249 :
1250 : ! tau_r and tau_g
1251 276 : IF (ASSOCIATED(tau_g_in)) THEN
1252 0 : nspins = SIZE(tau_g_in)
1253 0 : ALLOCATE (tau_g_out(nspins))
1254 0 : CALL qs_rho_set(rho_output, tau_g=tau_g_out)
1255 0 : IF (ASSOCIATED(tau_r_in)) THEN
1256 0 : ALLOCATE (tau_r_out(nspins))
1257 0 : CALL qs_rho_set(rho_output, tau_r=tau_r_out)
1258 : END IF
1259 0 : DO i = 1, nspins
1260 0 : CALL out_pw_pool%create_pw(tau_g_out(i))
1261 0 : CALL pw_transfer(tau_g_in(i), tau_g_out(i))
1262 0 : IF (ASSOCIATED(tau_r_in)) THEN
1263 0 : CALL out_pw_pool%create_pw(tau_r_out(i))
1264 0 : CALL pw_transfer(tau_g_out(i), tau_r_out(i))
1265 : END IF
1266 : END DO
1267 276 : ELSE IF (ASSOCIATED(tau_r_in)) THEN
1268 0 : nspins = SIZE(tau_r_in)
1269 0 : ALLOCATE (tau_r_out(nspins))
1270 0 : CALL qs_rho_set(rho_output, tau_r=tau_r_out)
1271 0 : DO i = 1, nspins
1272 0 : CALL out_pw_pool%create_pw(tau_r_out(i))
1273 0 : BLOCK
1274 : TYPE(pw_c1d_gs_type) :: gtau_in, gtau_out
1275 0 : CALL in_pw_pool%create_pw(gtau_in)
1276 0 : CALL out_pw_pool%create_pw(gtau_out)
1277 0 : CALL pw_transfer(tau_r_in(i), gtau_in)
1278 0 : CALL pw_transfer(gtau_in, gtau_out)
1279 0 : CALL pw_transfer(gtau_out, tau_r_out(i))
1280 0 : CALL out_pw_pool%give_back_pw(gtau_out)
1281 0 : CALL in_pw_pool%give_back_pw(gtau_in)
1282 : END BLOCK
1283 : END DO
1284 : END IF
1285 :
1286 : ! tot_rho_r
1287 276 : IF (ASSOCIATED(tot_rho_r_in)) THEN
1288 252 : nspins = SIZE(tot_rho_r_in)
1289 756 : ALLOCATE (tot_rho_r_out(nspins))
1290 252 : CALL qs_rho_set(rho_output, tot_rho_r=tot_rho_r_out)
1291 504 : DO i = 1, nspins
1292 504 : tot_rho_r_out(i) = tot_rho_r_in(i)
1293 : END DO
1294 : END IF
1295 :
1296 : ! tot_rho_g
1297 276 : IF (ASSOCIATED(tot_rho_g_in)) THEN
1298 0 : nspins = SIZE(tot_rho_g_in)
1299 0 : ALLOCATE (tot_rho_g_out(nspins))
1300 0 : CALL qs_rho_set(rho_output, tot_rho_g=tot_rho_g_out)
1301 0 : DO i = 1, nspins
1302 0 : tot_rho_g_out(i) = tot_rho_g_in(i)
1303 : END DO
1304 : END IF
1305 :
1306 : CALL qs_rho_set(rho_output, &
1307 : rho_g_valid=rho_g_valid_in, &
1308 : rho_r_valid=rho_r_valid_in, &
1309 : drho_g_valid=drho_g_valid_in, &
1310 : drho_r_valid=drho_r_valid_in, &
1311 : tau_r_valid=tau_r_valid_in, &
1312 : tau_g_valid=tau_g_valid_in, &
1313 : soft_valid=soft_valid_in, &
1314 276 : complex_rho_ao=complex_rho_ao)
1315 :
1316 276 : CALL timestop(handle)
1317 :
1318 276 : END SUBROUTINE qs_rho_transfer
1319 :
1320 : ! **************************************************************************************************
1321 : !> \brief rhoa(2) = alpha*rhoa(2)+beta*rhob(1)
1322 : !> \param rhoa ...
1323 : !> \param rhob ...
1324 : !> \param alpha ...
1325 : !> \param beta ...
1326 : ! **************************************************************************************************
1327 48 : SUBROUTINE qs_rho_scale_and_add_b(rhoa, rhob, alpha, beta)
1328 :
1329 : TYPE(qs_rho_type), INTENT(IN) :: rhoa, rhob
1330 : REAL(KIND=dp), INTENT(IN) :: alpha, beta
1331 :
1332 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_rho_scale_and_add_b'
1333 :
1334 : INTEGER :: handle
1335 48 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_g_a, tot_rho_g_b, tot_rho_r_a, &
1336 48 : tot_rho_r_b
1337 48 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_a, rho_ao_b, rho_ao_im_a, &
1338 48 : rho_ao_im_b
1339 48 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_a, rho_g_b, tau_g_a, tau_g_b
1340 48 : TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER :: drho_g_a, drho_g_b
1341 48 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_a, rho_r_b, tau_r_a, tau_r_b
1342 48 : TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER :: drho_r_a, drho_r_b
1343 : TYPE(pw_r3d_rs_type), POINTER :: rho_r_sccs_a, rho_r_sccs_b
1344 :
1345 48 : CALL timeset(routineN, handle)
1346 :
1347 48 : NULLIFY (rho_ao_a, rho_ao_im_a, rho_r_a, rho_g_a, drho_r_a, &
1348 48 : drho_g_a, tau_r_a, tau_g_a, tot_rho_r_a, tot_rho_g_a, rho_r_sccs_a)
1349 :
1350 : CALL qs_rho_get(rhoa, &
1351 : rho_ao=rho_ao_a, &
1352 : rho_ao_im=rho_ao_im_a, &
1353 : rho_r=rho_r_a, &
1354 : rho_g=rho_g_a, &
1355 : drho_r=drho_r_a, &
1356 : drho_g=drho_g_a, &
1357 : tau_r=tau_r_a, &
1358 : tau_g=tau_g_a, &
1359 : tot_rho_r=tot_rho_r_a, &
1360 : tot_rho_g=tot_rho_g_a, &
1361 48 : rho_r_sccs=rho_r_sccs_a)
1362 :
1363 48 : NULLIFY (rho_ao_b, rho_ao_im_b, rho_r_b, rho_g_b, drho_r_b, &
1364 48 : drho_g_b, tau_r_b, tau_g_b, tot_rho_r_b, tot_rho_g_b, rho_r_sccs_b)
1365 :
1366 : CALL qs_rho_get(rhob, &
1367 : rho_ao=rho_ao_b, &
1368 : rho_ao_im=rho_ao_im_b, &
1369 : rho_r=rho_r_b, &
1370 : rho_g=rho_g_b, &
1371 : drho_r=drho_r_b, &
1372 : drho_g=drho_g_b, &
1373 : tau_r=tau_r_b, &
1374 : tau_g=tau_g_b, &
1375 : tot_rho_r=tot_rho_r_b, &
1376 : tot_rho_g=tot_rho_g_b, &
1377 48 : rho_r_sccs=rho_r_sccs_b)
1378 : ! rho_ao
1379 48 : IF (ASSOCIATED(rho_ao_a) .AND. ASSOCIATED(rho_ao_b)) THEN
1380 48 : CALL dbcsr_add(rho_ao_a(2)%matrix, rho_ao_b(1)%matrix, alpha, beta)
1381 : END IF
1382 :
1383 : ! rho_ao_im
1384 48 : IF (ASSOCIATED(rho_ao_im_a) .AND. ASSOCIATED(rho_ao_im_b)) THEN
1385 0 : CALL dbcsr_add(rho_ao_im_a(2)%matrix, rho_ao_im_b(1)%matrix, alpha, beta)
1386 : END IF
1387 :
1388 : ! rho_r
1389 48 : IF (ASSOCIATED(rho_r_a) .AND. ASSOCIATED(rho_r_b)) THEN
1390 48 : CALL pw_axpy(rho_r_b(1), rho_r_a(2), beta, alpha)
1391 : END IF
1392 :
1393 : ! rho_g
1394 48 : IF (ASSOCIATED(rho_g_a) .AND. ASSOCIATED(rho_g_b)) THEN
1395 48 : CALL pw_axpy(rho_g_b(1), rho_g_a(2), beta, alpha)
1396 : END IF
1397 :
1398 : ! SCCS
1399 48 : IF (ASSOCIATED(rho_r_sccs_a) .AND. ASSOCIATED(rho_r_sccs_b)) THEN
1400 0 : CALL pw_axpy(rho_r_sccs_b, rho_r_sccs_a, beta, alpha)
1401 : END IF
1402 :
1403 : ! drho_r
1404 48 : IF (ASSOCIATED(drho_r_a) .AND. ASSOCIATED(drho_r_b)) THEN
1405 : CPASSERT(ASSOCIATED(drho_r_a) .AND. ASSOCIATED(drho_r_b)) ! not implemented
1406 : END IF
1407 :
1408 : ! drho_g
1409 48 : IF (ASSOCIATED(drho_g_a) .AND. ASSOCIATED(drho_g_b)) THEN
1410 : CPASSERT(ASSOCIATED(drho_g_a) .AND. ASSOCIATED(drho_g_b)) ! not implemented
1411 : END IF
1412 :
1413 : ! tau_r
1414 48 : IF (ASSOCIATED(tau_r_a) .AND. ASSOCIATED(tau_r_b)) THEN
1415 0 : CALL pw_axpy(tau_r_b(1), tau_r_a(2), beta, alpha)
1416 : END IF
1417 :
1418 : ! tau_g
1419 48 : IF (ASSOCIATED(tau_g_a) .AND. ASSOCIATED(tau_g_b)) THEN
1420 0 : CALL pw_axpy(tau_g_b(1), tau_g_a(2), beta, alpha)
1421 : END IF
1422 :
1423 : ! tot_rho_r
1424 48 : IF (ASSOCIATED(tot_rho_r_a) .AND. ASSOCIATED(tot_rho_r_b)) THEN
1425 0 : tot_rho_r_a(2) = alpha*tot_rho_r_a(2) + beta*tot_rho_r_b(1)
1426 : END IF
1427 :
1428 : ! tot_rho_g
1429 48 : IF (ASSOCIATED(tot_rho_g_a) .AND. ASSOCIATED(tot_rho_g_b)) THEN
1430 0 : tot_rho_g_a(2) = alpha*tot_rho_g_a(2) + beta*tot_rho_g_b(1)
1431 : END IF
1432 :
1433 48 : CALL timestop(handle)
1434 :
1435 48 : END SUBROUTINE qs_rho_scale_and_add_b
1436 :
1437 : ! **************************************************************************************************
1438 : !> \brief rhoa = alpha*rhoa+beta*rhob
1439 : !> \param rhoa ...
1440 : !> \param rhob ...
1441 : !> \param alpha ...
1442 : !> \param beta ...
1443 : ! **************************************************************************************************
1444 15696 : SUBROUTINE qs_rho_scale_and_add(rhoa, rhob, alpha, beta)
1445 :
1446 : TYPE(qs_rho_type), INTENT(IN) :: rhoa, rhob
1447 : REAL(KIND=dp), INTENT(IN) :: alpha, beta
1448 :
1449 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_rho_scale_and_add'
1450 :
1451 : INTEGER :: handle, i, j, nspina, nspinb, nspins
1452 15696 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_g_a, tot_rho_g_b, tot_rho_r_a, &
1453 15696 : tot_rho_r_b
1454 15696 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_a, rho_ao_b, rho_ao_im_a, &
1455 15696 : rho_ao_im_b
1456 15696 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_a, rho_g_b, tau_g_a, tau_g_b
1457 15696 : TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER :: drho_g_a, drho_g_b
1458 15696 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_a, rho_r_b, tau_r_a, tau_r_b
1459 15696 : TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER :: drho_r_a, drho_r_b
1460 : TYPE(pw_r3d_rs_type), POINTER :: rho_r_sccs_a, rho_r_sccs_b
1461 :
1462 15696 : CALL timeset(routineN, handle)
1463 :
1464 15696 : NULLIFY (rho_ao_a, rho_ao_im_a, rho_r_a, rho_g_a, drho_r_a, &
1465 15696 : drho_g_a, tau_r_a, tau_g_a, tot_rho_r_a, tot_rho_g_a, rho_r_sccs_a)
1466 :
1467 : CALL qs_rho_get(rhoa, &
1468 : rho_ao=rho_ao_a, &
1469 : rho_ao_im=rho_ao_im_a, &
1470 : rho_r=rho_r_a, &
1471 : rho_g=rho_g_a, &
1472 : drho_r=drho_r_a, &
1473 : drho_g=drho_g_a, &
1474 : tau_r=tau_r_a, &
1475 : tau_g=tau_g_a, &
1476 : tot_rho_r=tot_rho_r_a, &
1477 : tot_rho_g=tot_rho_g_a, &
1478 15696 : rho_r_sccs=rho_r_sccs_a)
1479 :
1480 15696 : NULLIFY (rho_ao_b, rho_ao_im_b, rho_r_b, rho_g_b, drho_r_b, &
1481 15696 : drho_g_b, tau_r_b, tau_g_b, tot_rho_r_b, tot_rho_g_b, rho_r_sccs_b)
1482 :
1483 : CALL qs_rho_get(rhob, &
1484 : rho_ao=rho_ao_b, &
1485 : rho_ao_im=rho_ao_im_b, &
1486 : rho_r=rho_r_b, &
1487 : rho_g=rho_g_b, &
1488 : drho_r=drho_r_b, &
1489 : drho_g=drho_g_b, &
1490 : tau_r=tau_r_b, &
1491 : tau_g=tau_g_b, &
1492 : tot_rho_r=tot_rho_r_b, &
1493 : tot_rho_g=tot_rho_g_b, &
1494 15696 : rho_r_sccs=rho_r_sccs_b)
1495 : ! rho_ao
1496 15696 : IF (ASSOCIATED(rho_ao_a) .AND. ASSOCIATED(rho_ao_b)) THEN
1497 15696 : nspina = SIZE(rho_ao_a)
1498 15696 : nspinb = SIZE(rho_ao_b)
1499 15696 : nspins = MIN(nspina, nspinb)
1500 33120 : DO i = 1, nspins
1501 33120 : CALL dbcsr_add(rho_ao_a(i)%matrix, rho_ao_b(i)%matrix, alpha, beta)
1502 : END DO
1503 : END IF
1504 :
1505 : ! rho_ao_im
1506 15696 : IF (ASSOCIATED(rho_ao_im_a) .AND. ASSOCIATED(rho_ao_im_b)) THEN
1507 0 : nspina = SIZE(rho_ao_im_a)
1508 0 : nspinb = SIZE(rho_ao_im_b)
1509 0 : nspins = MIN(nspina, nspinb)
1510 0 : DO i = 1, nspins
1511 0 : CALL dbcsr_add(rho_ao_im_a(i)%matrix, rho_ao_im_b(i)%matrix, alpha, beta)
1512 : END DO
1513 : END IF
1514 :
1515 : ! rho_r
1516 15696 : IF (ASSOCIATED(rho_r_a) .AND. ASSOCIATED(rho_r_b)) THEN
1517 15696 : nspina = SIZE(rho_ao_a)
1518 15696 : nspinb = SIZE(rho_ao_b)
1519 15696 : nspins = MIN(nspina, nspinb)
1520 33120 : DO i = 1, nspins
1521 33120 : CALL pw_axpy(rho_r_b(i), rho_r_a(i), beta, alpha)
1522 : END DO
1523 : END IF
1524 :
1525 : ! rho_g
1526 15696 : IF (ASSOCIATED(rho_g_a) .AND. ASSOCIATED(rho_g_b)) THEN
1527 15696 : nspina = SIZE(rho_ao_a)
1528 15696 : nspinb = SIZE(rho_ao_b)
1529 15696 : nspins = MIN(nspina, nspinb)
1530 33120 : DO i = 1, nspins
1531 33120 : CALL pw_axpy(rho_g_b(i), rho_g_a(i), beta, alpha)
1532 : END DO
1533 : END IF
1534 :
1535 : ! SCCS
1536 15696 : IF (ASSOCIATED(rho_r_sccs_a) .AND. ASSOCIATED(rho_r_sccs_b)) THEN
1537 0 : CALL pw_axpy(rho_r_sccs_b, rho_r_sccs_a, beta, alpha)
1538 : END IF
1539 :
1540 : ! drho_r
1541 15696 : IF (ASSOCIATED(drho_r_a) .AND. ASSOCIATED(drho_r_b)) THEN
1542 0 : CPASSERT(ALL(SHAPE(drho_r_a) == SHAPE(drho_r_b))) ! not implemented
1543 0 : DO j = 1, SIZE(drho_r_a, 2)
1544 0 : DO i = 1, SIZE(drho_r_a, 1)
1545 0 : CALL pw_axpy(drho_r_b(i, j), drho_r_a(i, j), beta, alpha)
1546 : END DO
1547 : END DO
1548 : END IF
1549 :
1550 : ! drho_g
1551 15696 : IF (ASSOCIATED(drho_g_a) .AND. ASSOCIATED(drho_g_b)) THEN
1552 0 : CPASSERT(ALL(SHAPE(drho_g_a) == SHAPE(drho_g_b))) ! not implemented
1553 0 : DO j = 1, SIZE(drho_g_a, 2)
1554 0 : DO i = 1, SIZE(drho_g_a, 1)
1555 0 : CALL pw_axpy(drho_g_b(i, j), drho_g_a(i, j), beta, alpha)
1556 : END DO
1557 : END DO
1558 : END IF
1559 :
1560 : ! tau_r
1561 15696 : IF (ASSOCIATED(tau_r_a) .AND. ASSOCIATED(tau_r_b)) THEN
1562 0 : nspina = SIZE(rho_ao_a)
1563 0 : nspinb = SIZE(rho_ao_b)
1564 0 : nspins = MIN(nspina, nspinb)
1565 0 : DO i = 1, nspins
1566 0 : CALL pw_axpy(tau_r_b(i), tau_r_a(i), beta, alpha)
1567 : END DO
1568 : END IF
1569 :
1570 : ! tau_g
1571 15696 : IF (ASSOCIATED(tau_g_a) .AND. ASSOCIATED(tau_g_b)) THEN
1572 0 : nspina = SIZE(rho_ao_a)
1573 0 : nspinb = SIZE(rho_ao_b)
1574 0 : nspins = MIN(nspina, nspinb)
1575 0 : DO i = 1, nspins
1576 0 : CALL pw_axpy(tau_g_b(i), tau_g_a(i), beta, alpha)
1577 : END DO
1578 : END IF
1579 :
1580 : ! tot_rho_r
1581 15696 : IF (ASSOCIATED(tot_rho_r_a) .AND. ASSOCIATED(tot_rho_r_b)) THEN
1582 0 : nspina = SIZE(rho_ao_a)
1583 0 : nspinb = SIZE(rho_ao_b)
1584 0 : nspins = MIN(nspina, nspinb)
1585 0 : DO i = 1, nspins
1586 0 : tot_rho_r_a(i) = alpha*tot_rho_r_a(i) + beta*tot_rho_r_b(i)
1587 : END DO
1588 : END IF
1589 :
1590 : ! tot_rho_g
1591 15696 : IF (ASSOCIATED(tot_rho_g_a) .AND. ASSOCIATED(tot_rho_g_b)) THEN
1592 0 : nspina = SIZE(rho_ao_a)
1593 0 : nspinb = SIZE(rho_ao_b)
1594 0 : nspins = MIN(nspina, nspinb)
1595 0 : DO i = 1, nspins
1596 0 : tot_rho_g_a(i) = alpha*tot_rho_g_a(i) + beta*tot_rho_g_b(i)
1597 : END DO
1598 : END IF
1599 :
1600 15696 : CALL timestop(handle)
1601 :
1602 15696 : END SUBROUTINE qs_rho_scale_and_add
1603 :
1604 : ! **************************************************************************************************
1605 : !> \brief Duplicates a pointer physically
1606 : !> \param rho_input The rho structure to be duplicated
1607 : !> \param rho_output The duplicate rho structure
1608 : !> \param qs_env The QS environment from which the auxiliary PW basis-set
1609 : !> pool is taken
1610 : !> \par History
1611 : !> 07.2005 initial create [tdk]
1612 : !> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch)
1613 : !> \note
1614 : !> Associated pointers are deallocated, nullified pointers are NOT accepted!
1615 : ! **************************************************************************************************
1616 8 : SUBROUTINE duplicate_rho_type(rho_input, rho_output, qs_env)
1617 :
1618 : TYPE(qs_rho_type), INTENT(INOUT) :: rho_input, rho_output
1619 : TYPE(qs_environment_type), POINTER :: qs_env
1620 :
1621 : CHARACTER(len=*), PARAMETER :: routineN = 'duplicate_rho_type'
1622 :
1623 : INTEGER :: handle, i, j, nspins
1624 : LOGICAL :: complex_rho_ao_in, drho_g_valid_in, drho_r_valid_in, rho_g_valid_in, &
1625 : rho_r_valid_in, soft_valid_in, tau_g_valid_in, tau_r_valid_in
1626 4 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_g_in, tot_rho_g_out, &
1627 4 : tot_rho_r_in, tot_rho_r_out
1628 4 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_im_in, rho_ao_im_out, rho_ao_in, &
1629 4 : rho_ao_out
1630 : TYPE(dft_control_type), POINTER :: dft_control
1631 4 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_in, rho_g_out, tau_g_in, tau_g_out
1632 4 : TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER :: drho_g_in, drho_g_out
1633 : TYPE(pw_env_type), POINTER :: pw_env
1634 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1635 4 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_in, rho_r_out, tau_r_in, tau_r_out
1636 4 : TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER :: drho_r_in, drho_r_out
1637 : TYPE(pw_r3d_rs_type), POINTER :: rho_r_sccs_in, rho_r_sccs_out
1638 :
1639 4 : CALL timeset(routineN, handle)
1640 :
1641 4 : NULLIFY (dft_control, pw_env, auxbas_pw_pool)
1642 4 : NULLIFY (rho_ao_in, rho_ao_out, rho_ao_im_in, rho_ao_im_out)
1643 4 : NULLIFY (rho_r_in, rho_r_out, rho_g_in, rho_g_out, drho_r_in, drho_r_out)
1644 4 : NULLIFY (drho_g_in, drho_g_out, tau_r_in, tau_r_out, tau_g_in, tau_g_out)
1645 4 : NULLIFY (tot_rho_r_in, tot_rho_r_out, tot_rho_g_in, tot_rho_g_out)
1646 4 : NULLIFY (rho_r_sccs_in, rho_r_sccs_out)
1647 :
1648 4 : CPASSERT(ASSOCIATED(qs_env))
1649 :
1650 4 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, dft_control=dft_control)
1651 4 : CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool)
1652 4 : nspins = dft_control%nspins
1653 :
1654 4 : CALL qs_rho_clear(rho_output)
1655 :
1656 : CALL qs_rho_get(rho_input, &
1657 : rho_ao=rho_ao_in, &
1658 : rho_ao_im=rho_ao_im_in, &
1659 : rho_r=rho_r_in, &
1660 : rho_g=rho_g_in, &
1661 : drho_r=drho_r_in, &
1662 : drho_g=drho_g_in, &
1663 : tau_r=tau_r_in, &
1664 : tau_g=tau_g_in, &
1665 : tot_rho_r=tot_rho_r_in, &
1666 : tot_rho_g=tot_rho_g_in, &
1667 : rho_g_valid=rho_g_valid_in, &
1668 : rho_r_valid=rho_r_valid_in, &
1669 : drho_g_valid=drho_g_valid_in, &
1670 : drho_r_valid=drho_r_valid_in, &
1671 : tau_r_valid=tau_r_valid_in, &
1672 : tau_g_valid=tau_g_valid_in, &
1673 : rho_r_sccs=rho_r_sccs_in, &
1674 : soft_valid=soft_valid_in, &
1675 4 : complex_rho_ao=complex_rho_ao_in)
1676 :
1677 : ! rho_ao
1678 4 : IF (ASSOCIATED(rho_ao_in)) THEN
1679 4 : CALL dbcsr_allocate_matrix_set(rho_ao_out, nspins)
1680 4 : CALL qs_rho_set(rho_output, rho_ao=rho_ao_out)
1681 8 : DO i = 1, nspins
1682 4 : ALLOCATE (rho_ao_out(i)%matrix)
1683 : CALL dbcsr_copy(rho_ao_out(i)%matrix, rho_ao_in(i)%matrix, &
1684 4 : name="myDensityMatrix_for_Spin_"//TRIM(ADJUSTL(cp_to_string(i))))
1685 8 : CALL dbcsr_set(rho_ao_out(i)%matrix, 0.0_dp)
1686 : END DO
1687 : END IF
1688 :
1689 : ! rho_ao_im
1690 4 : IF (ASSOCIATED(rho_ao_im_in)) THEN
1691 0 : CALL dbcsr_allocate_matrix_set(rho_ao_im_out, nspins)
1692 0 : CALL qs_rho_set(rho_output, rho_ao=rho_ao_im_out)
1693 0 : DO i = 1, nspins
1694 0 : ALLOCATE (rho_ao_im_out(i)%matrix)
1695 : CALL dbcsr_copy(rho_ao_im_out(i)%matrix, rho_ao_im_in(i)%matrix, &
1696 0 : name="myImagDensityMatrix_for_Spin_"//TRIM(ADJUSTL(cp_to_string(i))))
1697 0 : CALL dbcsr_set(rho_ao_im_out(i)%matrix, 0.0_dp)
1698 : END DO
1699 : END IF
1700 :
1701 : ! rho_r
1702 4 : IF (ASSOCIATED(rho_r_in)) THEN
1703 16 : ALLOCATE (rho_r_out(nspins))
1704 4 : CALL qs_rho_set(rho_output, rho_r=rho_r_out)
1705 8 : DO i = 1, nspins
1706 4 : CALL auxbas_pw_pool%create_pw(rho_r_out(i))
1707 8 : CALL pw_copy(rho_r_in(i), rho_r_out(i))
1708 : END DO
1709 : END IF
1710 :
1711 : ! rho_g
1712 4 : IF (ASSOCIATED(rho_g_in)) THEN
1713 16 : ALLOCATE (rho_g_out(nspins))
1714 4 : CALL qs_rho_set(rho_output, rho_g=rho_g_out)
1715 8 : DO i = 1, nspins
1716 4 : CALL auxbas_pw_pool%create_pw(rho_g_out(i))
1717 8 : CALL pw_copy(rho_g_in(i), rho_g_out(i))
1718 : END DO
1719 : END IF
1720 :
1721 : ! SCCS
1722 4 : IF (ASSOCIATED(rho_r_sccs_in)) THEN
1723 0 : CALL qs_rho_set(rho_output, rho_r_sccs=rho_r_sccs_out)
1724 0 : CALL auxbas_pw_pool%create_pw(rho_r_sccs_out)
1725 0 : CALL pw_copy(rho_r_sccs_in, rho_r_sccs_out)
1726 : END IF
1727 :
1728 : ! drho_r and drho_g are only needed if calculated by collocation
1729 4 : IF (dft_control%drho_by_collocation) THEN
1730 : ! drho_r
1731 0 : IF (ASSOCIATED(drho_r_in)) THEN
1732 0 : ALLOCATE (drho_r_out(3, nspins))
1733 0 : CALL qs_rho_set(rho_output, drho_r=drho_r_out)
1734 0 : DO j = 1, nspins
1735 0 : DO i = 1, 3
1736 0 : CALL auxbas_pw_pool%create_pw(drho_r_out(i, j))
1737 0 : CALL pw_copy(drho_r_in(i, j), drho_r_out(i, j))
1738 : END DO
1739 : END DO
1740 : END IF
1741 :
1742 : ! drho_g
1743 0 : IF (ASSOCIATED(drho_g_in)) THEN
1744 0 : ALLOCATE (drho_g_out(3, nspins))
1745 0 : CALL qs_rho_set(rho_output, drho_g=drho_g_out)
1746 0 : DO j = 1, nspins
1747 0 : DO i = 1, 3
1748 0 : CALL auxbas_pw_pool%create_pw(drho_g_out(i, j))
1749 0 : CALL pw_copy(drho_g_in(i, j), drho_g_out(i, j))
1750 : END DO
1751 : END DO
1752 : END IF
1753 : END IF
1754 :
1755 : ! tau_r and tau_g are only needed in the case of Meta-GGA XC-functionals
1756 : ! are used. Therefore they are only allocated if
1757 : ! dft_control%use_kinetic_energy_density is true
1758 4 : IF (dft_control%use_kinetic_energy_density) THEN
1759 : ! tau_r
1760 0 : IF (ASSOCIATED(tau_r_in)) THEN
1761 0 : ALLOCATE (tau_r_out(nspins))
1762 0 : CALL qs_rho_set(rho_output, tau_r=tau_r_out)
1763 0 : DO i = 1, nspins
1764 0 : CALL auxbas_pw_pool%create_pw(tau_r_out(i))
1765 0 : CALL pw_copy(tau_r_in(i), tau_r_out(i))
1766 : END DO
1767 : END IF
1768 :
1769 : ! tau_g
1770 0 : IF (ASSOCIATED(tau_g_in)) THEN
1771 0 : ALLOCATE (tau_g_out(nspins))
1772 0 : CALL qs_rho_set(rho_output, tau_g=tau_g_out)
1773 0 : DO i = 1, nspins
1774 0 : CALL auxbas_pw_pool%create_pw(tau_g_out(i))
1775 0 : CALL pw_copy(tau_g_in(i), tau_g_out(i))
1776 : END DO
1777 : END IF
1778 : END IF
1779 :
1780 : CALL qs_rho_set(rho_output, &
1781 : rho_g_valid=rho_g_valid_in, &
1782 : rho_r_valid=rho_r_valid_in, &
1783 : drho_g_valid=drho_g_valid_in, &
1784 : drho_r_valid=drho_r_valid_in, &
1785 : tau_r_valid=tau_r_valid_in, &
1786 : tau_g_valid=tau_g_valid_in, &
1787 : soft_valid=soft_valid_in, &
1788 4 : complex_rho_ao=complex_rho_ao_in)
1789 :
1790 : ! tot_rho_r
1791 4 : IF (ASSOCIATED(tot_rho_r_in)) THEN
1792 12 : ALLOCATE (tot_rho_r_out(nspins))
1793 4 : CALL qs_rho_set(rho_output, tot_rho_r=tot_rho_r_out)
1794 8 : DO i = 1, nspins
1795 8 : tot_rho_r_out(i) = tot_rho_r_in(i)
1796 : END DO
1797 : END IF
1798 :
1799 : ! tot_rho_g
1800 4 : IF (ASSOCIATED(tot_rho_g_in)) THEN
1801 0 : ALLOCATE (tot_rho_g_out(nspins))
1802 0 : CALL qs_rho_set(rho_output, tot_rho_g=tot_rho_g_out)
1803 0 : DO i = 1, nspins
1804 0 : tot_rho_g_out(i) = tot_rho_g_in(i)
1805 : END DO
1806 :
1807 : END IF
1808 :
1809 4 : CALL timestop(handle)
1810 :
1811 4 : END SUBROUTINE duplicate_rho_type
1812 :
1813 : ! **************************************************************************************************
1814 : !> \brief (Re-)allocates rho_ao_im from real part rho_ao
1815 : !> \param rho ...
1816 : !> \param qs_env ...
1817 : ! **************************************************************************************************
1818 116 : SUBROUTINE allocate_rho_ao_imag_from_real(rho, qs_env)
1819 : TYPE(qs_rho_type), POINTER :: rho
1820 : TYPE(qs_environment_type), POINTER :: qs_env
1821 :
1822 : CHARACTER(LEN=default_string_length) :: headline
1823 : INTEGER :: i, ic, nimages, nspins
1824 116 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_im_kp, rho_ao_kp
1825 : TYPE(dbcsr_type), POINTER :: template
1826 : TYPE(dft_control_type), POINTER :: dft_control
1827 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1828 116 : POINTER :: sab_orb
1829 :
1830 116 : NULLIFY (rho_ao_im_kp, rho_ao_kp, dft_control, template, sab_orb)
1831 :
1832 : CALL get_qs_env(qs_env, &
1833 : dft_control=dft_control, &
1834 116 : sab_orb=sab_orb)
1835 :
1836 116 : CALL qs_rho_get(rho, rho_ao_im_kp=rho_ao_im_kp, rho_ao_kp=rho_ao_kp)
1837 :
1838 116 : nspins = dft_control%nspins
1839 116 : nimages = dft_control%nimages
1840 :
1841 116 : CPASSERT(nspins == SIZE(rho_ao_kp, 1))
1842 116 : CPASSERT(nimages == SIZE(rho_ao_kp, 2))
1843 :
1844 116 : CALL dbcsr_allocate_matrix_set(rho_ao_im_kp, nspins, nimages)
1845 116 : CALL qs_rho_set(rho, rho_ao_im_kp=rho_ao_im_kp)
1846 254 : DO i = 1, nspins
1847 392 : DO ic = 1, nimages
1848 138 : IF (nspins > 1) THEN
1849 44 : IF (i == 1) THEN
1850 22 : headline = "IMAGINARY PART OF DENSITY MATRIX FOR ALPHA SPIN"
1851 : ELSE
1852 22 : headline = "IMAGINARY PART OF DENSITY MATRIX FOR BETA SPIN"
1853 : END IF
1854 : ELSE
1855 94 : headline = "IMAGINARY PART OF DENSITY MATRIX"
1856 : END IF
1857 138 : ALLOCATE (rho_ao_im_kp(i, ic)%matrix)
1858 138 : template => rho_ao_kp(i, ic)%matrix ! base on real part, but anti-symmetric
1859 : CALL dbcsr_create(matrix=rho_ao_im_kp(i, ic)%matrix, template=template, &
1860 138 : name=TRIM(headline), matrix_type=dbcsr_type_antisymmetric)
1861 138 : CALL cp_dbcsr_alloc_block_from_nbl(rho_ao_im_kp(i, ic)%matrix, sab_orb)
1862 276 : CALL dbcsr_set(rho_ao_im_kp(i, ic)%matrix, 0.0_dp)
1863 : END DO
1864 : END DO
1865 :
1866 116 : END SUBROUTINE allocate_rho_ao_imag_from_real
1867 :
1868 : END MODULE qs_rho_methods
|
