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