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