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