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