Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief contains the structure
10 : !> \par History
11 : !> 11.2003 created [fawzi]
12 : !> \author fawzi
13 : ! **************************************************************************************************
14 : MODULE xc_rho_set_types
15 : USE cp_array_utils, ONLY: cp_3d_r_cp_type
16 : USE kinds, ONLY: dp
17 : USE pw_grid_types, ONLY: pw_grid_type
18 : USE pw_methods, ONLY: pw_copy, &
19 : pw_transfer
20 : USE pw_pool_types, ONLY: &
21 : pw_pool_type
22 : USE pw_spline_utils, ONLY: pw_spline_scale_deriv
23 : USE pw_types, ONLY: &
24 : pw_c1d_gs_type, &
25 : pw_r3d_rs_type
26 : USE xc_input_constants, ONLY: xc_deriv_pw, &
27 : xc_deriv_spline2, &
28 : xc_deriv_spline3, &
29 : xc_rho_no_smooth
30 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_equal, &
31 : xc_rho_cflags_setall, &
32 : xc_rho_cflags_type
33 : USE xc_util, ONLY: xc_pw_gradient, &
34 : xc_pw_laplace, &
35 : xc_pw_smooth
36 : #include "../base/base_uses.f90"
37 :
38 : IMPLICIT NONE
39 : PRIVATE
40 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
41 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_rho_set_types'
42 :
43 : PUBLIC :: xc_rho_set_type
44 : PUBLIC :: xc_rho_set_create, xc_rho_set_release, &
45 : xc_rho_set_update, xc_rho_set_get, xc_rho_set_recover_pw
46 :
47 : ! **************************************************************************************************
48 : !> \brief represent a density, with all the representation and data needed
49 : !> to perform a functional evaluation
50 : !> \param local_bounds the part of the 3d array on which the functional is
51 : !> computed
52 : !> \param owns which components are owned by this structure (and should
53 : !> be deallocated
54 : !> \param has which components are present and up to date
55 : !> \param rho the density
56 : !> \param drho the gradient of the density (x,y and z direction)
57 : !> \param norm_drho the norm of the gradient of the density
58 : !> \param rhoa , rhob: spin alpha and beta parts of the density in the LSD case
59 : !> \param drhoa , drhob: gradient of the spin alpha and beta parts of the density
60 : !> in the LSD case (x,y and z direction)
61 : !> \param norm_drhoa , norm_drhob: norm of the gradient of rhoa and rhob
62 : !> \param rho_ 1_3: rho^(1.0_dp/3.0_dp)
63 : !> \param rhoa_ 1_3, rhob_1_3: rhoa^(1.0_dp/3.0_dp), rhob^(1.0_dp/3.0_dp)
64 : !> \param tau the kinetic (KohnSham) part of rho
65 : !> \param tau_a the kinetic (KohnSham) part of rhoa
66 : !> \param tau_b the kinetic (KohnSham) part of rhob
67 : !> \note
68 : !> the use of 3d arrays is the result of trying to use only basic
69 : !> types (to be generic and independent from the method), and
70 : !> avoiding copies using the actual structure.
71 : !> only the part defined by local bounds is guaranteed to be present,
72 : !> and it is the only meaningful part.
73 : !> \par History
74 : !> 11.2003 created [fawzi & thomas]
75 : !> 12.2008 added laplace parts [mguidon]
76 : !> \author fawzi & thomas
77 : ! **************************************************************************************************
78 : TYPE xc_rho_set_type
79 : INTEGER, DIMENSION(2, 3) :: local_bounds = -1
80 : REAL(kind=dp) :: rho_cutoff = EPSILON(0.0_dp), drho_cutoff = EPSILON(0.0_dp), tau_cutoff = EPSILON(0.0_dp)
81 : TYPE(xc_rho_cflags_type) :: owns = xc_rho_cflags_type(), has = xc_rho_cflags_type()
82 : ! for spin restricted systems
83 : REAL(KIND=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: rho => NULL()
84 : TYPE(cp_3d_r_cp_type), DIMENSION(3) :: drho = cp_3d_r_cp_type()
85 : REAL(KIND=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: norm_drho => NULL()
86 : REAL(KIND=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: rho_1_3 => NULL()
87 : REAL(kind=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: tau => NULL()
88 : ! for UNrestricted systems
89 : REAL(KIND=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: rhoa => NULL(), rhob => NULL()
90 : TYPE(cp_3d_r_cp_type), DIMENSION(3) :: drhoa = cp_3d_r_cp_type(), &
91 : drhob = cp_3d_r_cp_type()
92 : REAL(KIND=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: norm_drhoa => NULL(), norm_drhob => NULL()
93 : REAL(kind=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: rhoa_1_3 => NULL(), rhob_1_3 => NULL()
94 : REAL(kind=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: tau_a => NULL(), tau_b => NULL()
95 : REAL(kind=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: laplace_rho => NULL(), laplace_rhoa => NULL(), &
96 : laplace_rhob => NULL()
97 : END TYPE xc_rho_set_type
98 :
99 : CONTAINS
100 :
101 : ! **************************************************************************************************
102 : !> \brief allocates and does (minimal) initialization of a rho_set
103 : !> \param rho_set the structure to allocate
104 : !> \param local_bounds ...
105 : !> \param rho_cutoff ...
106 : !> \param drho_cutoff ...
107 : !> \param tau_cutoff ...
108 : ! **************************************************************************************************
109 269097 : SUBROUTINE xc_rho_set_create(rho_set, local_bounds, rho_cutoff, drho_cutoff, &
110 : tau_cutoff)
111 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
112 : INTEGER, DIMENSION(2, 3), INTENT(in) :: local_bounds
113 : REAL(kind=dp), INTENT(in), OPTIONAL :: rho_cutoff, drho_cutoff, tau_cutoff
114 :
115 269097 : IF (PRESENT(rho_cutoff)) rho_set%rho_cutoff = rho_cutoff
116 269097 : IF (PRESENT(drho_cutoff)) rho_set%drho_cutoff = drho_cutoff
117 269097 : IF (PRESENT(tau_cutoff)) rho_set%tau_cutoff = tau_cutoff
118 2690970 : rho_set%local_bounds = local_bounds
119 269097 : CALL xc_rho_cflags_setall(rho_set%owns, .TRUE.)
120 269097 : CALL xc_rho_cflags_setall(rho_set%has, .FALSE.)
121 269097 : END SUBROUTINE xc_rho_set_create
122 :
123 : ! **************************************************************************************************
124 : !> \brief releases the given rho_set
125 : !> \param rho_set the structure to release
126 : !> \param pw_pool the plae where to give back the arrays
127 : ! **************************************************************************************************
128 269097 : SUBROUTINE xc_rho_set_release(rho_set, pw_pool)
129 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
130 : TYPE(pw_pool_type), OPTIONAL, POINTER :: pw_pool
131 :
132 : INTEGER :: i
133 :
134 269097 : IF (PRESENT(pw_pool)) THEN
135 126427 : IF (ASSOCIATED(pw_pool)) THEN
136 126427 : CALL xc_rho_set_clean(rho_set, pw_pool)
137 : ELSE
138 0 : CPABORT("pw_pool must be associated")
139 : END IF
140 : END IF
141 :
142 1076388 : rho_set%local_bounds(1, :) = -HUGE(0) ! we want to crash...
143 1076388 : rho_set%local_bounds(1, :) = HUGE(0)
144 269097 : IF (rho_set%owns%rho .AND. ASSOCIATED(rho_set%rho)) THEN
145 123673 : DEALLOCATE (rho_set%rho)
146 : ELSE
147 145424 : NULLIFY (rho_set%rho)
148 : END IF
149 269097 : IF (rho_set%owns%rho_spin) THEN
150 120410 : IF (ASSOCIATED(rho_set%rhoa)) THEN
151 18907 : DEALLOCATE (rho_set%rhoa)
152 : END IF
153 120410 : IF (ASSOCIATED(rho_set%rhob)) THEN
154 18803 : DEALLOCATE (rho_set%rhob)
155 : END IF
156 : ELSE
157 148687 : NULLIFY (rho_set%rhoa, rho_set%rhob)
158 : END IF
159 269097 : IF (rho_set%owns%rho_1_3 .AND. ASSOCIATED(rho_set%rho_1_3)) THEN
160 5281 : DEALLOCATE (rho_set%rho_1_3)
161 : ELSE
162 263816 : NULLIFY (rho_set%rho_1_3)
163 : END IF
164 269097 : IF (rho_set%owns%rho_spin) THEN
165 120410 : IF (ASSOCIATED(rho_set%rhoa_1_3)) THEN
166 2452 : DEALLOCATE (rho_set%rhoa_1_3)
167 : END IF
168 120410 : IF (ASSOCIATED(rho_set%rhob_1_3)) THEN
169 2452 : DEALLOCATE (rho_set%rhob_1_3)
170 : END IF
171 : ELSE
172 148687 : NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
173 : END IF
174 269097 : IF (rho_set%owns%drho) THEN
175 516584 : DO i = 1, 3
176 516584 : IF (ASSOCIATED(rho_set%drho(i)%array)) THEN
177 190962 : DEALLOCATE (rho_set%drho(i)%array)
178 : END IF
179 : END DO
180 : ELSE
181 559804 : DO i = 1, 3
182 559804 : NULLIFY (rho_set%drho(i)%array)
183 : END DO
184 : END IF
185 269097 : IF (rho_set%owns%drho_spin) THEN
186 473920 : DO i = 1, 3
187 355440 : IF (ASSOCIATED(rho_set%drhoa(i)%array)) THEN
188 32316 : DEALLOCATE (rho_set%drhoa(i)%array)
189 : END IF
190 473920 : IF (ASSOCIATED(rho_set%drhob(i)%array)) THEN
191 32160 : DEALLOCATE (rho_set%drhob(i)%array)
192 : END IF
193 : END DO
194 : ELSE
195 602468 : DO i = 1, 3
196 602468 : NULLIFY (rho_set%drhoa(i)%array, rho_set%drhob(i)%array)
197 : END DO
198 : END IF
199 269097 : IF (rho_set%owns%laplace_rho .AND. ASSOCIATED(rho_set%laplace_rho)) THEN
200 202 : DEALLOCATE (rho_set%laplace_rho)
201 : ELSE
202 268895 : NULLIFY (rho_set%laplace_rho)
203 : END IF
204 :
205 269097 : IF (rho_set%owns%norm_drho .AND. ASSOCIATED(rho_set%norm_drho)) THEN
206 96528 : DEALLOCATE (rho_set%norm_drho)
207 : ELSE
208 172569 : NULLIFY (rho_set%norm_drho)
209 : END IF
210 269097 : IF (rho_set%owns%laplace_rho_spin) THEN
211 116690 : IF (ASSOCIATED(rho_set%laplace_rhoa)) THEN
212 50 : DEALLOCATE (rho_set%laplace_rhoa)
213 : END IF
214 116690 : IF (ASSOCIATED(rho_set%laplace_rhob)) THEN
215 50 : DEALLOCATE (rho_set%laplace_rhob)
216 : END IF
217 : ELSE
218 152407 : NULLIFY (rho_set%laplace_rhoa, rho_set%laplace_rhob)
219 : END IF
220 :
221 269097 : IF (rho_set%owns%norm_drho_spin) THEN
222 118480 : IF (ASSOCIATED(rho_set%norm_drhoa)) THEN
223 11651 : DEALLOCATE (rho_set%norm_drhoa)
224 : END IF
225 118480 : IF (ASSOCIATED(rho_set%norm_drhob)) THEN
226 11599 : DEALLOCATE (rho_set%norm_drhob)
227 : END IF
228 : ELSE
229 150617 : NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
230 : END IF
231 269097 : IF (rho_set%owns%tau .AND. ASSOCIATED(rho_set%tau)) THEN
232 1092 : DEALLOCATE (rho_set%tau)
233 : ELSE
234 268005 : NULLIFY (rho_set%tau)
235 : END IF
236 269097 : IF (rho_set%owns%tau_spin) THEN
237 116640 : IF (ASSOCIATED(rho_set%tau_a)) THEN
238 2 : DEALLOCATE (rho_set%tau_a)
239 : END IF
240 116640 : IF (ASSOCIATED(rho_set%tau_b)) THEN
241 2 : DEALLOCATE (rho_set%tau_b)
242 : END IF
243 : ELSE
244 152457 : NULLIFY (rho_set%tau_a, rho_set%tau_b)
245 : END IF
246 269097 : END SUBROUTINE xc_rho_set_release
247 :
248 : ! **************************************************************************************************
249 : !> \brief returns the various attributes of rho_set
250 : !> \param rho_set the object you want info about
251 : !> \param can_return_null if true the object returned can be null,
252 : !> if false (the default) it stops with an error if a requested
253 : !> component is not associated
254 : !> \param rho ...
255 : !> \param drho ...
256 : !> \param norm_drho ...
257 : !> \param rhoa ...
258 : !> \param rhob ...
259 : !> \param norm_drhoa ...
260 : !> \param norm_drhob ...
261 : !> \param rho_1_3 ...
262 : !> \param rhoa_1_3 ...
263 : !> \param rhob_1_3 ...
264 : !> \param laplace_rho ...
265 : !> \param laplace_rhoa ...
266 : !> \param laplace_rhob ...
267 : !> \param drhoa ...
268 : !> \param drhob ...
269 : !> \param rho_cutoff ...
270 : !> \param drho_cutoff ...
271 : !> \param tau_cutoff ...
272 : !> \param tau ...
273 : !> \param tau_a ...
274 : !> \param tau_b ...
275 : !> \param local_bounds ...
276 : ! **************************************************************************************************
277 1071804 : SUBROUTINE xc_rho_set_get(rho_set, can_return_null, rho, drho, norm_drho, &
278 : rhoa, rhob, norm_drhoa, norm_drhob, rho_1_3, rhoa_1_3, &
279 : rhob_1_3, laplace_rho, laplace_rhoa, laplace_rhob, drhoa, drhob, rho_cutoff, &
280 : drho_cutoff, tau_cutoff, tau, tau_a, tau_b, local_bounds)
281 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
282 : LOGICAL, INTENT(in), OPTIONAL :: can_return_null
283 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
284 : POINTER :: rho
285 : TYPE(cp_3d_r_cp_type), DIMENSION(3), OPTIONAL :: drho
286 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
287 : POINTER :: norm_drho, rhoa, rhob, norm_drhoa, &
288 : norm_drhob, rho_1_3, rhoa_1_3, &
289 : rhob_1_3, laplace_rho, laplace_rhoa, &
290 : laplace_rhob
291 : TYPE(cp_3d_r_cp_type), DIMENSION(3), OPTIONAL :: drhoa, drhob
292 : REAL(kind=dp), INTENT(out), OPTIONAL :: rho_cutoff, drho_cutoff, tau_cutoff
293 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
294 : POINTER :: tau, tau_a, tau_b
295 : INTEGER, DIMENSION(2, 3), INTENT(OUT), OPTIONAL :: local_bounds
296 :
297 : INTEGER :: i
298 : LOGICAL :: my_can_return_null
299 :
300 1071804 : my_can_return_null = .FALSE.
301 1071804 : IF (PRESENT(can_return_null)) my_can_return_null = can_return_null
302 :
303 1071804 : IF (PRESENT(rho)) THEN
304 277832 : rho => rho_set%rho
305 277832 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rho))
306 : END IF
307 1071804 : IF (PRESENT(drho)) THEN
308 737428 : DO i = 1, 3
309 553071 : drho(i)%array => rho_set%drho(i)%array
310 737428 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rho_set%drho(i)%array))
311 : END DO
312 : END IF
313 1071804 : IF (PRESENT(norm_drho)) THEN
314 388930 : norm_drho => rho_set%norm_drho
315 388930 : CPASSERT(my_can_return_null .OR. ASSOCIATED(norm_drho))
316 : END IF
317 1071804 : IF (PRESENT(laplace_rho)) THEN
318 13320 : laplace_rho => rho_set%laplace_rho
319 13320 : CPASSERT(my_can_return_null .OR. ASSOCIATED(laplace_rho))
320 : END IF
321 1071804 : IF (PRESENT(rhoa)) THEN
322 150874 : rhoa => rho_set%rhoa
323 150874 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rhoa))
324 : END IF
325 1071804 : IF (PRESENT(rhob)) THEN
326 150460 : rhob => rho_set%rhob
327 150460 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rhob))
328 : END IF
329 1071804 : IF (PRESENT(drhoa)) THEN
330 613764 : DO i = 1, 3
331 460323 : drhoa(i)%array => rho_set%drhoa(i)%array
332 613764 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rho_set%drhoa(i)%array))
333 : END DO
334 : END IF
335 1071804 : IF (PRESENT(drhob)) THEN
336 612716 : DO i = 1, 3
337 459537 : drhob(i)%array => rho_set%drhob(i)%array
338 612716 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rho_set%drhob(i)%array))
339 : END DO
340 : END IF
341 1071804 : IF (PRESENT(laplace_rhoa)) THEN
342 2570 : laplace_rhoa => rho_set%laplace_rhoa
343 2570 : CPASSERT(my_can_return_null .OR. ASSOCIATED(laplace_rhoa))
344 : END IF
345 1071804 : IF (PRESENT(laplace_rhob)) THEN
346 2570 : laplace_rhob => rho_set%laplace_rhob
347 2570 : CPASSERT(my_can_return_null .OR. ASSOCIATED(laplace_rhob))
348 : END IF
349 1071804 : IF (PRESENT(norm_drhoa)) THEN
350 205061 : norm_drhoa => rho_set%norm_drhoa
351 205061 : CPASSERT(my_can_return_null .OR. ASSOCIATED(norm_drhoa))
352 : END IF
353 1071804 : IF (PRESENT(norm_drhob)) THEN
354 205061 : norm_drhob => rho_set%norm_drhob
355 205061 : CPASSERT(my_can_return_null .OR. ASSOCIATED(norm_drhob))
356 : END IF
357 1071804 : IF (PRESENT(rho_1_3)) THEN
358 21277 : rho_1_3 => rho_set%rho_1_3
359 21277 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rho_1_3))
360 : END IF
361 1071804 : IF (PRESENT(rhoa_1_3)) THEN
362 3348 : rhoa_1_3 => rho_set%rhoa_1_3
363 3348 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rhoa_1_3))
364 : END IF
365 1071804 : IF (PRESENT(rhob_1_3)) THEN
366 3348 : rhob_1_3 => rho_set%rhob_1_3
367 3348 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rhob_1_3))
368 : END IF
369 1071804 : IF (PRESENT(tau)) THEN
370 15092 : tau => rho_set%tau
371 15092 : CPASSERT(my_can_return_null .OR. ASSOCIATED(tau))
372 : END IF
373 1071804 : IF (PRESENT(tau_a)) THEN
374 2754 : tau_a => rho_set%tau_a
375 2754 : CPASSERT(my_can_return_null .OR. ASSOCIATED(tau_a))
376 : END IF
377 1071804 : IF (PRESENT(tau_b)) THEN
378 2754 : tau_b => rho_set%tau_b
379 2754 : CPASSERT(my_can_return_null .OR. ASSOCIATED(tau_b))
380 : END IF
381 1071804 : IF (PRESENT(rho_cutoff)) rho_cutoff = rho_set%rho_cutoff
382 1071804 : IF (PRESENT(drho_cutoff)) drho_cutoff = rho_set%drho_cutoff
383 1071804 : IF (PRESENT(tau_cutoff)) tau_cutoff = rho_set%tau_cutoff
384 2753544 : IF (PRESENT(local_bounds)) local_bounds = rho_set%local_bounds
385 1071804 : END SUBROUTINE xc_rho_set_get
386 :
387 : #:mute
388 : #:def recover(variable)
389 : #! Determine component of the actual data
390 : #:set var_long_pw = (variable+"(i)" if variable.startswith("drho") else variable)
391 : #:set var_long_rho = (variable+"(i)%array" if variable.startswith("drho") else variable)
392 : #! Determine the flag name
393 : #! Remove spin states and potential underscore
394 : #:set is_1_3 = variable.endswith("_1_3")
395 : #:set var_base = variable.strip("_13")
396 : #:set is_spin = var_base.endswith("a") or var_base.endswith("b")
397 : #:set var_base = var_base.strip("ab_")
398 : #:set var_cflags = (var_base if not is_spin else var_base+"_spin")
399 : #:set var_cflags = (var_cflags if not is_1_3 else var_cflags+"_1_3")
400 : IF (PRESENT(${variable}$)) THEN
401 : #:if variable.startswith("drho")
402 : DO i = 1, 3
403 : #:else
404 : NULLIFY (${var_long_pw}$)
405 : ALLOCATE (${var_long_pw}$)
406 : #:endif
407 : CALL xc_rho_set_recover_pw_low(${var_long_pw}$, rho_set%${var_long_rho}$, pw_grid, pw_pool#{if variable =="drho"}#, rho_set%drhoa(i)%array, rho_set%drhob(i)%array#{endif}#)
408 : #:if not variable.startswith("drho")
409 : NULLIFY (rho_set%${var_long_rho}$)
410 : #:else
411 : END DO
412 : #:endif
413 : owns_data = #{if variable =="drho"}#.TRUE.#{else}#rho_set%owns%${var_cflags}$#{endif}#
414 : END IF
415 : #:enddef
416 : #:endmute
417 :
418 : ! **************************************************************************************************
419 : !> \brief Shifts association of the requested array to a pw grid
420 : !> Requires that the corresponding component of rho_set is associated
421 : !> If owns_data returns TRUE, the caller has to allocate the data later
422 : !> It is allowed to task for only one component per call.
423 : !> In case of drho, the array is allocated if not internally available and calculated from drhoa and drhob.
424 : !> \param rho_set the object you want info about
425 : !> \param pw_grid ...
426 : !> \param pw_pool ...
427 : !> \param owns_data ...
428 : !> \param rho ...
429 : !> \param drho ...
430 : !> \param norm_drho ...
431 : !> \param rhoa ...
432 : !> \param rhob ...
433 : !> \param norm_drhoa ...
434 : !> \param norm_drhob ...
435 : !> \param rho_1_3 ...
436 : !> \param rhoa_1_3 ...
437 : !> \param rhob_1_3 ...
438 : !> \param laplace_rho ...
439 : !> \param laplace_rhoa ...
440 : !> \param laplace_rhob ...
441 : !> \param drhoa ...
442 : !> \param drhob ...
443 : !> \param tau ...
444 : !> \param tau_a ...
445 : !> \param tau_b ...
446 : ! **************************************************************************************************
447 528756 : SUBROUTINE xc_rho_set_recover_pw(rho_set, pw_grid, pw_pool, owns_data, rho, drho, norm_drho, &
448 : rhoa, rhob, norm_drhoa, norm_drhob, rho_1_3, rhoa_1_3, &
449 : rhob_1_3, laplace_rho, laplace_rhoa, laplace_rhob, drhoa, drhob, &
450 : tau, tau_a, tau_b)
451 : TYPE(xc_rho_set_type) :: rho_set
452 : TYPE(pw_r3d_rs_type), DIMENSION(3), OPTIONAL, INTENT(OUT) :: drho, drhoa, drhob
453 : TYPE(pw_r3d_rs_type), OPTIONAL, POINTER :: rho, norm_drho, rhoa, rhob, norm_drhoa, &
454 : norm_drhob, rho_1_3, rhoa_1_3, &
455 : rhob_1_3, laplace_rho, laplace_rhoa, &
456 : laplace_rhob, tau, tau_a, tau_b
457 : TYPE(pw_grid_type), POINTER, INTENT(IN) :: pw_grid
458 : TYPE(pw_pool_type), POINTER, INTENT(IN) :: pw_pool
459 : LOGICAL, INTENT(OUT) :: owns_data
460 :
461 : INTEGER :: i
462 :
463 : #:for variable in ["rho", "drho", "norm_drho", "rhoa", "rhob", "norm_drhoa", "norm_drhob", "rho_1_3", "rhoa_1_3", "rhob_1_3", "laplace_rho", "laplace_rhoa", "laplace_rhob", "drhoa", "drhob", "tau", "tau_a", "tau_b"]
464 452795 : $:recover(variable)
465 : #:endfor
466 :
467 90559 : END SUBROUTINE xc_rho_set_recover_pw
468 :
469 271677 : SUBROUTINE xc_rho_set_recover_pw_low(rho_pw, rho, pw_grid, pw_pool, rhoa, rhob)
470 : TYPE(pw_r3d_rs_type), INTENT(OUT) :: rho_pw
471 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER, CONTIGUOUS :: rho
472 : TYPE(pw_grid_type), POINTER, INTENT(IN) :: pw_grid
473 : TYPE(pw_pool_type), POINTER, INTENT(IN) :: pw_pool
474 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER, OPTIONAL :: rhoa, rhob
475 :
476 271677 : IF (ASSOCIATED(rho)) THEN
477 230157 : CALL rho_pw%create(pw_grid=pw_grid, array_ptr=rho)
478 230157 : NULLIFY (rho)
479 41520 : ELSE IF (PRESENT(rhoa) .AND. PRESENT(rhob)) THEN
480 41520 : CALL pw_pool%create_pw(rho_pw)
481 41520 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(rho_pw,rhoa,rhob)
482 : rho_pw%array(:, :, :) = rhoa(:, :, :) + rhob(:, :, :)
483 : !$OMP END PARALLEL WORKSHARE
484 : ELSE
485 : CALL cp_abort(__LOCATION__, "Either component or its spin parts (if applicable) "// &
486 0 : "have to be associated in rho_set!")
487 : END IF
488 :
489 271677 : END SUBROUTINE xc_rho_set_recover_pw_low
490 :
491 : ! **************************************************************************************************
492 : !> \brief cleans (releases) most of the data stored in the rho_set giving back
493 : !> what it can to the pw_pool
494 : !> \param rho_set the rho_set to be cleaned
495 : !> \param pw_pool place to give back 3d arrays,...
496 : !> \author Fawzi Mohamed
497 : ! **************************************************************************************************
498 286128 : SUBROUTINE xc_rho_set_clean(rho_set, pw_pool)
499 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
500 : TYPE(pw_pool_type), POINTER :: pw_pool
501 :
502 : INTEGER :: idir
503 :
504 286128 : IF (rho_set%owns%rho) THEN
505 256657 : CALL pw_pool%give_back_cr3d(rho_set%rho)
506 : ELSE
507 29471 : NULLIFY (rho_set%rho)
508 : END IF
509 286128 : IF (rho_set%owns%rho_1_3) THEN
510 161173 : CALL pw_pool%give_back_cr3d(rho_set%rho_1_3)
511 : ELSE
512 124955 : NULLIFY (rho_set%rho_1_3)
513 : END IF
514 286128 : IF (rho_set%owns%drho) THEN
515 820728 : DO idir = 1, 3
516 820728 : CALL pw_pool%give_back_cr3d(rho_set%drho(idir)%array)
517 : END DO
518 : ELSE
519 323784 : DO idir = 1, 3
520 323784 : NULLIFY (rho_set%drho(idir)%array)
521 : END DO
522 : END IF
523 286128 : IF (rho_set%owns%norm_drho) THEN
524 225850 : CALL pw_pool%give_back_cr3d(rho_set%norm_drho)
525 : ELSE
526 60278 : NULLIFY (rho_set%norm_drho)
527 : END IF
528 286128 : IF (rho_set%owns%laplace_rho) THEN
529 153233 : CALL pw_pool%give_back_cr3d(rho_set%laplace_rho)
530 : ELSE
531 132895 : NULLIFY (rho_set%laplace_rho)
532 : END IF
533 286128 : IF (rho_set%owns%tau) THEN
534 152457 : CALL pw_pool%give_back_cr3d(rho_set%tau)
535 : ELSE
536 133671 : NULLIFY (rho_set%tau)
537 : END IF
538 286128 : IF (rho_set%owns%rho_spin) THEN
539 181928 : CALL pw_pool%give_back_cr3d(rho_set%rhoa)
540 181928 : CALL pw_pool%give_back_cr3d(rho_set%rhob)
541 : ELSE
542 104200 : NULLIFY (rho_set%rhoa, rho_set%rhob)
543 : END IF
544 286128 : IF (rho_set%owns%rho_spin_1_3) THEN
545 154201 : CALL pw_pool%give_back_cr3d(rho_set%rhoa_1_3)
546 154201 : CALL pw_pool%give_back_cr3d(rho_set%rhob_1_3)
547 : ELSE
548 131927 : NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
549 : END IF
550 286128 : IF (rho_set%owns%drho_spin) THEN
551 672492 : DO idir = 1, 3
552 504369 : CALL pw_pool%give_back_cr3d(rho_set%drhoa(idir)%array)
553 672492 : CALL pw_pool%give_back_cr3d(rho_set%drhob(idir)%array)
554 : END DO
555 : ELSE
556 472020 : DO idir = 1, 3
557 472020 : NULLIFY (rho_set%drhoa(idir)%array, rho_set%drhob(idir)%array)
558 : END DO
559 : END IF
560 286128 : IF (rho_set%owns%laplace_rho_spin) THEN
561 152757 : CALL pw_pool%give_back_cr3d(rho_set%laplace_rhoa)
562 152757 : CALL pw_pool%give_back_cr3d(rho_set%laplace_rhob)
563 : ELSE
564 133371 : NULLIFY (rho_set%laplace_rhoa, rho_set%laplace_rhob)
565 : END IF
566 286128 : IF (rho_set%owns%norm_drho_spin) THEN
567 171003 : CALL pw_pool%give_back_cr3d(rho_set%norm_drhoa)
568 171003 : CALL pw_pool%give_back_cr3d(rho_set%norm_drhob)
569 : ELSE
570 115125 : NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
571 : END IF
572 286128 : IF (rho_set%owns%tau_spin) THEN
573 152457 : CALL pw_pool%give_back_cr3d(rho_set%tau_a)
574 152457 : CALL pw_pool%give_back_cr3d(rho_set%tau_b)
575 : ELSE
576 133671 : NULLIFY (rho_set%tau_a, rho_set%tau_b)
577 : END IF
578 :
579 286128 : CALL xc_rho_cflags_setall(rho_set%has, .FALSE.)
580 286128 : CALL xc_rho_cflags_setall(rho_set%owns, .FALSE.)
581 :
582 286128 : END SUBROUTINE xc_rho_set_clean
583 :
584 : ! **************************************************************************************************
585 : !> \brief updates the given rho set with the density given by
586 : !> rho_r (and rho_g). The rho set will contain the components specified
587 : !> in needs
588 : !> \param rho_set the rho_set to update
589 : !> \param rho_r the new density (in r space)
590 : !> \param rho_g the new density (in g space, needed for some
591 : !> derivatives)
592 : !> \param tau ...
593 : !> \param needs the components of rho that are needed
594 : !> \param xc_deriv_method_id ...
595 : !> \param xc_rho_smooth_id ...
596 : !> \param pw_pool pool for the allocation of pw and array
597 : ! **************************************************************************************************
598 159701 : SUBROUTINE xc_rho_set_update(rho_set, rho_r, rho_g, tau, needs, &
599 : xc_deriv_method_id, xc_rho_smooth_id, pw_pool, spinflip)
600 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
601 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: rho_r
602 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
603 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER, INTENT(IN) :: tau
604 : TYPE(xc_rho_cflags_type), INTENT(in) :: needs
605 : INTEGER, INTENT(IN) :: xc_deriv_method_id, xc_rho_smooth_id
606 : TYPE(pw_pool_type), POINTER :: pw_pool
607 : LOGICAL, OPTIONAL :: spinflip
608 :
609 : REAL(KIND=dp), PARAMETER :: f13 = (1.0_dp/3.0_dp)
610 :
611 : INTEGER :: i, idir, ispin, j, k, nspins
612 : LOGICAL :: gradient_f, my_rho_g_local, &
613 : needs_laplace, needs_rho_g, do_sf
614 : REAL(kind=dp) :: rho_cutoff
615 479103 : TYPE(pw_r3d_rs_type), DIMENSION(2) :: laplace_rho_r
616 1437309 : TYPE(pw_r3d_rs_type), DIMENSION(3, 2) :: drho_r
617 : TYPE(pw_c1d_gs_type) :: my_rho_g, tmp_g
618 479103 : TYPE(pw_r3d_rs_type), DIMENSION(2) :: my_rho_r
619 :
620 159701 : do_sf = .FALSE.
621 14120 : IF (PRESENT(spinflip)) do_sf = spinflip
622 :
623 1597010 : IF (ANY(rho_set%local_bounds /= pw_pool%pw_grid%bounds_local)) &
624 0 : CPABORT("pw_pool cr3d have different size than expected")
625 159701 : nspins = SIZE(rho_r)
626 1597010 : rho_set%local_bounds = rho_r(1)%pw_grid%bounds_local
627 159701 : rho_cutoff = 0.5*rho_set%rho_cutoff
628 :
629 159701 : my_rho_g_local = .FALSE.
630 : ! some checks
631 126564 : SELECT CASE (nspins)
632 : CASE (1)
633 126564 : IF (.NOT. do_sf) THEN
634 126460 : CPASSERT(.NOT. needs%rho_spin)
635 126460 : CPASSERT(.NOT. needs%drho_spin)
636 126460 : CPASSERT(.NOT. needs%norm_drho_spin)
637 126460 : CPASSERT(.NOT. needs%rho_spin_1_3)
638 126460 : CPASSERT(.NOT. needs%tau_spin)
639 126460 : CPASSERT(.NOT. needs%laplace_rho_spin)
640 : ELSE
641 104 : CPASSERT(.NOT. needs%rho)
642 104 : CPASSERT(.NOT. needs%drho)
643 104 : CPASSERT(.NOT. needs%rho_1_3)
644 104 : CPASSERT(.NOT. needs%tau)
645 104 : CPASSERT(.NOT. needs%laplace_rho)
646 : END IF
647 : CASE (2)
648 33137 : CPASSERT(.NOT. needs%rho)
649 33137 : CPASSERT(.NOT. needs%drho)
650 33137 : CPASSERT(.NOT. needs%rho_1_3)
651 33137 : CPASSERT(.NOT. needs%tau)
652 33137 : CPASSERT(.NOT. needs%laplace_rho)
653 : CASE default
654 159701 : CPABORT("Unknown number of spin states")
655 : END SELECT
656 :
657 159701 : CALL xc_rho_set_clean(rho_set, pw_pool=pw_pool)
658 :
659 159701 : needs_laplace = (needs%laplace_rho .OR. needs%laplace_rho_spin)
660 : gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
661 : needs%drho .OR. needs%norm_drho .OR. &
662 159701 : needs_laplace)
663 : needs_rho_g = ((xc_deriv_method_id == xc_deriv_spline3 .OR. &
664 : xc_deriv_method_id == xc_deriv_spline2 .OR. &
665 159701 : xc_deriv_method_id == xc_deriv_pw)) .AND. (gradient_f .OR. needs_laplace)
666 159701 : IF ((gradient_f .AND. needs_laplace) .AND. &
667 : (xc_deriv_method_id /= xc_deriv_pw)) THEN
668 : CALL cp_abort(__LOCATION__, &
669 : "MGGA functionals that require the Laplacian are "// &
670 0 : "only compatible with 'XC_DERIV PW' and 'XC_SMOOTH_RHO NONE'")
671 : END IF
672 :
673 159701 : IF (needs_rho_g) THEN
674 87263 : CALL pw_pool%create_pw(tmp_g)
675 : END IF
676 352539 : DO ispin = 1, nspins
677 192838 : CALL pw_pool%create_pw(my_rho_r(ispin))
678 : ! introduce a smoothing kernel on the density
679 192838 : IF (xc_rho_smooth_id == xc_rho_no_smooth) THEN
680 192298 : IF (needs_rho_g) THEN
681 107057 : IF (ASSOCIATED(rho_g)) THEN
682 91141 : my_rho_g_local = .FALSE.
683 91141 : my_rho_g = rho_g(ispin)
684 : END IF
685 : END IF
686 :
687 192298 : CALL pw_copy(rho_r(ispin), my_rho_r(ispin))
688 : ELSE
689 540 : CALL xc_pw_smooth(rho_r(ispin), my_rho_r(ispin), xc_rho_smooth_id)
690 : END IF
691 :
692 352539 : IF (gradient_f) THEN ! calculate the grad of rho
693 : ! normally when you need the gradient you need the whole gradient
694 : ! (for the partial integration)
695 : ! deriv rho
696 436948 : DO idir = 1, 3
697 436948 : CALL pw_pool%create_pw(drho_r(idir, ispin))
698 : END DO
699 109237 : IF (needs_rho_g) THEN
700 107135 : IF (.NOT. ASSOCIATED(my_rho_g%pw_grid)) THEN
701 15994 : my_rho_g_local = .TRUE.
702 15994 : CALL pw_pool%create_pw(my_rho_g)
703 15994 : CALL pw_transfer(my_rho_r(ispin), my_rho_g)
704 : END IF
705 91141 : IF (.NOT. my_rho_g_local .AND. (xc_deriv_method_id == xc_deriv_spline2 .OR. &
706 : xc_deriv_method_id == xc_deriv_spline3)) THEN
707 7122 : CALL pw_pool%create_pw(my_rho_g)
708 7122 : my_rho_g_local = .TRUE.
709 7122 : CALL pw_copy(rho_g(ispin), my_rho_g)
710 : END IF
711 : END IF
712 109237 : IF (needs%laplace_rho .OR. needs%laplace_rho_spin) THEN
713 1678 : CALL pw_pool%create_pw(laplace_rho_r(ispin))
714 1678 : CALL xc_pw_laplace(my_rho_g, pw_pool, xc_deriv_method_id, laplace_rho_r(ispin), tmp_g=tmp_g)
715 : END IF
716 109237 : CALL xc_pw_gradient(my_rho_r(ispin), my_rho_g, tmp_g, drho_r(:, ispin), xc_deriv_method_id)
717 :
718 109237 : IF (needs_rho_g) THEN
719 107135 : IF (my_rho_g_local) THEN
720 23116 : my_rho_g_local = .FALSE.
721 23116 : CALL pw_pool%give_back_pw(my_rho_g)
722 : END IF
723 : END IF
724 :
725 109237 : IF (xc_deriv_method_id /= xc_deriv_pw) THEN
726 9290 : CALL pw_spline_scale_deriv(drho_r(:, ispin))
727 : END IF
728 :
729 : END IF
730 :
731 : END DO
732 :
733 159701 : IF (ASSOCIATED(tmp_g%pw_grid)) THEN
734 87263 : CALL pw_pool%give_back_pw(tmp_g)
735 : END IF
736 :
737 126564 : SELECT CASE (nspins)
738 : CASE (1)
739 126564 : IF (.NOT. do_sf) THEN
740 126460 : IF (needs%rho_1_3) THEN
741 9942 : CALL pw_pool%create_cr3d(rho_set%rho_1_3)
742 9942 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,my_rho_r)
743 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
744 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
745 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
746 : rho_set%rho_1_3(i, j, k) = MAX(my_rho_r(1)%array(i, j, k), 0.0_dp)**f13
747 : END DO
748 : END DO
749 : END DO
750 9942 : rho_set%owns%rho_1_3 = .TRUE.
751 9942 : rho_set%has%rho_1_3 = .TRUE.
752 : END IF
753 126460 : IF (needs%rho) THEN
754 126460 : rho_set%rho => my_rho_r(1)%array
755 126460 : NULLIFY (my_rho_r(1)%array)
756 126460 : rho_set%owns%rho = .TRUE.
757 126460 : rho_set%has%rho = .TRUE.
758 : END IF
759 126460 : IF (needs%norm_drho) THEN
760 68517 : CALL pw_pool%create_cr3d(rho_set%norm_drho)
761 68517 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
762 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
763 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
764 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
765 : rho_set%norm_drho(i, j, k) = SQRT( &
766 : drho_r(1, 1)%array(i, j, k)**2 + &
767 : drho_r(2, 1)%array(i, j, k)**2 + &
768 : drho_r(3, 1)%array(i, j, k)**2)
769 : END DO
770 : END DO
771 : END DO
772 68517 : rho_set%owns%norm_drho = .TRUE.
773 68517 : rho_set%has%norm_drho = .TRUE.
774 : END IF
775 126460 : IF (needs%laplace_rho) THEN
776 978 : rho_set%laplace_rho => laplace_rho_r(1)%array
777 978 : NULLIFY (laplace_rho_r(1)%array)
778 978 : rho_set%owns%laplace_rho = .TRUE.
779 978 : rho_set%has%laplace_rho = .TRUE.
780 : END IF
781 :
782 126460 : IF (needs%drho) THEN
783 260924 : DO idir = 1, 3
784 195693 : rho_set%drho(idir)%array => drho_r(idir, 1)%array
785 260924 : NULLIFY (drho_r(idir, 1)%array)
786 : END DO
787 65231 : rho_set%owns%drho = .TRUE.
788 65231 : rho_set%has%drho = .TRUE.
789 : END IF
790 : ELSE
791 104 : IF (needs%norm_drho) THEN
792 52 : CALL pw_pool%create_cr3d(rho_set%norm_drho)
793 52 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
794 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
795 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
796 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
797 : rho_set%norm_drho(i, j, k) = SQRT( &
798 : drho_r(1, 1)%array(i, j, k)**2 + &
799 : drho_r(2, 1)%array(i, j, k)**2 + &
800 : drho_r(3, 1)%array(i, j, k)**2)
801 : END DO
802 : END DO
803 : END DO
804 52 : rho_set%owns%norm_drho = .TRUE.
805 52 : rho_set%has%norm_drho = .TRUE.
806 : END IF
807 104 : IF (needs%rho_spin) THEN
808 :
809 104 : rho_set%rhoa => my_rho_r(1)%array
810 104 : NULLIFY (my_rho_r(1)%array)
811 :
812 104 : rho_set%owns%rho_spin = .TRUE.
813 104 : rho_set%has%rho_spin = .TRUE.
814 : END IF
815 104 : IF (needs%norm_drho_spin) THEN
816 52 : CALL pw_pool%create_cr3d(rho_set%norm_drhoa)
817 52 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
818 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
819 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
820 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
821 : rho_set%norm_drhoa(i, j, k) = SQRT( &
822 : drho_r(1, 1)%array(i, j, k)**2 + &
823 : drho_r(2, 1)%array(i, j, k)**2 + &
824 : drho_r(3, 1)%array(i, j, k)**2)
825 : END DO
826 : END DO
827 : END DO
828 52 : rho_set%owns%norm_drho_spin = .TRUE.
829 52 : rho_set%has%norm_drho_spin = .TRUE.
830 : END IF
831 104 : IF (needs%laplace_rho_spin) THEN
832 0 : rho_set%laplace_rhoa => laplace_rho_r(1)%array
833 0 : NULLIFY (laplace_rho_r(1)%array)
834 :
835 0 : rho_set%owns%laplace_rho_spin = .TRUE.
836 0 : rho_set%has%laplace_rho_spin = .TRUE.
837 : END IF
838 104 : IF (needs%drho_spin) THEN
839 208 : DO idir = 1, 3
840 156 : rho_set%drhoa(idir)%array => drho_r(idir, 1)%array
841 208 : NULLIFY (drho_r(idir, 1)%array)
842 : END DO
843 52 : rho_set%owns%drho_spin = .TRUE.
844 52 : rho_set%has%drho_spin = .TRUE.
845 : END IF
846 : END IF
847 : CASE (2)
848 33137 : IF (needs%rho_spin_1_3) THEN
849 1756 : CALL pw_pool%create_cr3d(rho_set%rhoa_1_3)
850 : !assume that the bounds are the same?
851 1756 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,my_rho_r)
852 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
853 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
854 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
855 : rho_set%rhoa_1_3(i, j, k) = MAX(my_rho_r(1)%array(i, j, k), 0.0_dp)**f13
856 : END DO
857 : END DO
858 : END DO
859 1756 : CALL pw_pool%create_cr3d(rho_set%rhob_1_3)
860 : !assume that the bounds are the same?
861 1756 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,my_rho_r)
862 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
863 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
864 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
865 : rho_set%rhob_1_3(i, j, k) = MAX(my_rho_r(2)%array(i, j, k), 0.0_dp)**f13
866 : END DO
867 : END DO
868 : END DO
869 1756 : rho_set%owns%rho_spin_1_3 = .TRUE.
870 1756 : rho_set%has%rho_spin_1_3 = .TRUE.
871 : END IF
872 33137 : IF (needs%norm_drho) THEN
873 :
874 19150 : CALL pw_pool%create_cr3d(rho_set%norm_drho)
875 19150 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
876 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
877 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
878 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
879 : rho_set%norm_drho(i, j, k) = SQRT( &
880 : (drho_r(1, 1)%array(i, j, k) + drho_r(1, 2)%array(i, j, k))**2 + &
881 : (drho_r(2, 1)%array(i, j, k) + drho_r(2, 2)%array(i, j, k))**2 + &
882 : (drho_r(3, 1)%array(i, j, k) + drho_r(3, 2)%array(i, j, k))**2)
883 : END DO
884 : END DO
885 : END DO
886 :
887 19150 : rho_set%owns%norm_drho = .TRUE.
888 19150 : rho_set%has%norm_drho = .TRUE.
889 : END IF
890 33137 : IF (needs%rho_spin) THEN
891 :
892 33137 : rho_set%rhoa => my_rho_r(1)%array
893 33137 : NULLIFY (my_rho_r(1)%array)
894 :
895 33137 : rho_set%rhob => my_rho_r(2)%array
896 33137 : NULLIFY (my_rho_r(2)%array)
897 :
898 33137 : rho_set%owns%rho_spin = .TRUE.
899 33137 : rho_set%has%rho_spin = .TRUE.
900 : END IF
901 33137 : IF (needs%norm_drho_spin) THEN
902 :
903 20334 : CALL pw_pool%create_cr3d(rho_set%norm_drhoa)
904 20334 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
905 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
906 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
907 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
908 : rho_set%norm_drhoa(i, j, k) = SQRT( &
909 : drho_r(1, 1)%array(i, j, k)**2 + &
910 : drho_r(2, 1)%array(i, j, k)**2 + &
911 : drho_r(3, 1)%array(i, j, k)**2)
912 : END DO
913 : END DO
914 : END DO
915 :
916 20334 : CALL pw_pool%create_cr3d(rho_set%norm_drhob)
917 20334 : rho_set%owns%norm_drho_spin = .TRUE.
918 20334 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
919 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
920 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
921 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
922 : rho_set%norm_drhob(i, j, k) = SQRT( &
923 : drho_r(1, 2)%array(i, j, k)**2 + &
924 : drho_r(2, 2)%array(i, j, k)**2 + &
925 : drho_r(3, 2)%array(i, j, k)**2)
926 : END DO
927 : END DO
928 : END DO
929 :
930 20334 : rho_set%owns%norm_drho_spin = .TRUE.
931 20334 : rho_set%has%norm_drho_spin = .TRUE.
932 : END IF
933 33137 : IF (needs%laplace_rho_spin) THEN
934 350 : rho_set%laplace_rhoa => laplace_rho_r(1)%array
935 350 : NULLIFY (laplace_rho_r(1)%array)
936 :
937 350 : rho_set%laplace_rhob => laplace_rho_r(2)%array
938 350 : NULLIFY (laplace_rho_r(2)%array)
939 :
940 350 : rho_set%owns%laplace_rho_spin = .TRUE.
941 350 : rho_set%has%laplace_rho_spin = .TRUE.
942 : END IF
943 192838 : IF (needs%drho_spin) THEN
944 69816 : DO idir = 1, 3
945 52362 : rho_set%drhoa(idir)%array => drho_r(idir, 1)%array
946 52362 : NULLIFY (drho_r(idir, 1)%array)
947 52362 : rho_set%drhob(idir)%array => drho_r(idir, 2)%array
948 69816 : NULLIFY (drho_r(idir, 2)%array)
949 : END DO
950 17454 : rho_set%owns%drho_spin = .TRUE.
951 17454 : rho_set%has%drho_spin = .TRUE.
952 : END IF
953 : END SELECT
954 : ! post cleanup
955 352539 : DO ispin = 1, nspins
956 192838 : IF (needs%laplace_rho .OR. needs%laplace_rho_spin) THEN
957 1678 : CALL pw_pool%give_back_pw(laplace_rho_r(ispin))
958 : END IF
959 931053 : DO idir = 1, 3
960 771352 : CALL pw_pool%give_back_pw(drho_r(idir, ispin))
961 : END DO
962 : END DO
963 352539 : DO ispin = 1, nspins
964 352539 : CALL pw_pool%give_back_pw(my_rho_r(ispin))
965 : END DO
966 :
967 : ! tau part
968 159701 : IF (needs%tau .OR. needs%tau_spin) THEN
969 3790 : CPASSERT(ASSOCIATED(tau))
970 8392 : DO ispin = 1, nspins
971 164303 : CPASSERT(ASSOCIATED(tau(ispin)%array))
972 : END DO
973 : END IF
974 159701 : IF (needs%tau) THEN
975 2978 : rho_set%tau => tau(1)%array
976 2978 : rho_set%owns%tau = .FALSE.
977 2978 : rho_set%has%tau = .TRUE.
978 : END IF
979 159701 : IF (needs%tau_spin) THEN
980 812 : rho_set%tau_a => tau(1)%array
981 812 : rho_set%tau_b => tau(2)%array
982 812 : rho_set%owns%tau_spin = .FALSE.
983 812 : rho_set%has%tau_spin = .TRUE.
984 : END IF
985 :
986 159701 : CPASSERT(xc_rho_cflags_equal(rho_set%has, needs))
987 :
988 319402 : END SUBROUTINE xc_rho_set_update
989 :
990 0 : END MODULE xc_rho_set_types
|
