(git:1145456)
Loading...
Searching...
No Matches
xc_rho_set_types.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief contains the structure
10!> \par History
11!> 11.2003 created [fawzi]
12!> \author fawzi
13! **************************************************************************************************
14MODULE 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
99CONTAINS
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 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 IF (PRESENT(rho_cutoff)) rho_set%rho_cutoff = rho_cutoff
116 IF (PRESENT(drho_cutoff)) rho_set%drho_cutoff = drho_cutoff
117 IF (PRESENT(tau_cutoff)) rho_set%tau_cutoff = tau_cutoff
118 rho_set%local_bounds = local_bounds
119 CALL xc_rho_cflags_setall(rho_set%owns, .true.)
120 CALL xc_rho_cflags_setall(rho_set%has, .false.)
121 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 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 IF (PRESENT(pw_pool)) THEN
135 IF (ASSOCIATED(pw_pool)) THEN
136 CALL xc_rho_set_clean(rho_set, pw_pool)
137 ELSE
138 cpabort("pw_pool must be associated")
139 END IF
140 END IF
141
142 rho_set%local_bounds(1, :) = -huge(0) ! we want to crash...
143 rho_set%local_bounds(1, :) = huge(0)
144 IF (rho_set%owns%rho .AND. ASSOCIATED(rho_set%rho)) THEN
145 DEALLOCATE (rho_set%rho)
146 ELSE
147 NULLIFY (rho_set%rho)
148 END IF
149 IF (rho_set%owns%rho_spin) THEN
150 IF (ASSOCIATED(rho_set%rhoa)) THEN
151 DEALLOCATE (rho_set%rhoa)
152 END IF
153 IF (ASSOCIATED(rho_set%rhob)) THEN
154 DEALLOCATE (rho_set%rhob)
155 END IF
156 ELSE
157 NULLIFY (rho_set%rhoa, rho_set%rhob)
158 END IF
159 IF (rho_set%owns%rho_1_3 .AND. ASSOCIATED(rho_set%rho_1_3)) THEN
160 DEALLOCATE (rho_set%rho_1_3)
161 ELSE
162 NULLIFY (rho_set%rho_1_3)
163 END IF
164 IF (rho_set%owns%rho_spin) THEN
165 IF (ASSOCIATED(rho_set%rhoa_1_3)) THEN
166 DEALLOCATE (rho_set%rhoa_1_3)
167 END IF
168 IF (ASSOCIATED(rho_set%rhob_1_3)) THEN
169 DEALLOCATE (rho_set%rhob_1_3)
170 END IF
171 ELSE
172 NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
173 END IF
174 IF (rho_set%owns%drho) THEN
175 DO i = 1, 3
176 IF (ASSOCIATED(rho_set%drho(i)%array)) THEN
177 DEALLOCATE (rho_set%drho(i)%array)
178 END IF
179 END DO
180 ELSE
181 DO i = 1, 3
182 NULLIFY (rho_set%drho(i)%array)
183 END DO
184 END IF
185 IF (rho_set%owns%drho_spin) THEN
186 DO i = 1, 3
187 IF (ASSOCIATED(rho_set%drhoa(i)%array)) THEN
188 DEALLOCATE (rho_set%drhoa(i)%array)
189 END IF
190 IF (ASSOCIATED(rho_set%drhob(i)%array)) THEN
191 DEALLOCATE (rho_set%drhob(i)%array)
192 END IF
193 END DO
194 ELSE
195 DO i = 1, 3
196 NULLIFY (rho_set%drhoa(i)%array, rho_set%drhob(i)%array)
197 END DO
198 END IF
199 IF (rho_set%owns%laplace_rho .AND. ASSOCIATED(rho_set%laplace_rho)) THEN
200 DEALLOCATE (rho_set%laplace_rho)
201 ELSE
202 NULLIFY (rho_set%laplace_rho)
203 END IF
204
205 IF (rho_set%owns%norm_drho .AND. ASSOCIATED(rho_set%norm_drho)) THEN
206 DEALLOCATE (rho_set%norm_drho)
207 ELSE
208 NULLIFY (rho_set%norm_drho)
209 END IF
210 IF (rho_set%owns%laplace_rho_spin) THEN
211 IF (ASSOCIATED(rho_set%laplace_rhoa)) THEN
212 DEALLOCATE (rho_set%laplace_rhoa)
213 END IF
214 IF (ASSOCIATED(rho_set%laplace_rhob)) THEN
215 DEALLOCATE (rho_set%laplace_rhob)
216 END IF
217 ELSE
218 NULLIFY (rho_set%laplace_rhoa, rho_set%laplace_rhob)
219 END IF
220
221 IF (rho_set%owns%norm_drho_spin) THEN
222 IF (ASSOCIATED(rho_set%norm_drhoa)) THEN
223 DEALLOCATE (rho_set%norm_drhoa)
224 END IF
225 IF (ASSOCIATED(rho_set%norm_drhob)) THEN
226 DEALLOCATE (rho_set%norm_drhob)
227 END IF
228 ELSE
229 NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
230 END IF
231 IF (rho_set%owns%tau .AND. ASSOCIATED(rho_set%tau)) THEN
232 DEALLOCATE (rho_set%tau)
233 ELSE
234 NULLIFY (rho_set%tau)
235 END IF
236 IF (rho_set%owns%tau_spin) THEN
237 IF (ASSOCIATED(rho_set%tau_a)) THEN
238 DEALLOCATE (rho_set%tau_a)
239 END IF
240 IF (ASSOCIATED(rho_set%tau_b)) THEN
241 DEALLOCATE (rho_set%tau_b)
242 END IF
243 ELSE
244 NULLIFY (rho_set%tau_a, rho_set%tau_b)
245 END IF
246 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 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 my_can_return_null = .false.
301 IF (PRESENT(can_return_null)) my_can_return_null = can_return_null
302
303 IF (PRESENT(rho)) THEN
304 rho => rho_set%rho
305 cpassert(my_can_return_null .OR. ASSOCIATED(rho))
306 END IF
307 IF (PRESENT(drho)) THEN
308 DO i = 1, 3
309 drho(i)%array => rho_set%drho(i)%array
310 cpassert(my_can_return_null .OR. ASSOCIATED(rho_set%drho(i)%array))
311 END DO
312 END IF
313 IF (PRESENT(norm_drho)) THEN
314 norm_drho => rho_set%norm_drho
315 cpassert(my_can_return_null .OR. ASSOCIATED(norm_drho))
316 END IF
317 IF (PRESENT(laplace_rho)) THEN
318 laplace_rho => rho_set%laplace_rho
319 cpassert(my_can_return_null .OR. ASSOCIATED(laplace_rho))
320 END IF
321 IF (PRESENT(rhoa)) THEN
322 rhoa => rho_set%rhoa
323 cpassert(my_can_return_null .OR. ASSOCIATED(rhoa))
324 END IF
325 IF (PRESENT(rhob)) THEN
326 rhob => rho_set%rhob
327 cpassert(my_can_return_null .OR. ASSOCIATED(rhob))
328 END IF
329 IF (PRESENT(drhoa)) THEN
330 DO i = 1, 3
331 drhoa(i)%array => rho_set%drhoa(i)%array
332 cpassert(my_can_return_null .OR. ASSOCIATED(rho_set%drhoa(i)%array))
333 END DO
334 END IF
335 IF (PRESENT(drhob)) THEN
336 DO i = 1, 3
337 drhob(i)%array => rho_set%drhob(i)%array
338 cpassert(my_can_return_null .OR. ASSOCIATED(rho_set%drhob(i)%array))
339 END DO
340 END IF
341 IF (PRESENT(laplace_rhoa)) THEN
342 laplace_rhoa => rho_set%laplace_rhoa
343 cpassert(my_can_return_null .OR. ASSOCIATED(laplace_rhoa))
344 END IF
345 IF (PRESENT(laplace_rhob)) THEN
346 laplace_rhob => rho_set%laplace_rhob
347 cpassert(my_can_return_null .OR. ASSOCIATED(laplace_rhob))
348 END IF
349 IF (PRESENT(norm_drhoa)) THEN
350 norm_drhoa => rho_set%norm_drhoa
351 cpassert(my_can_return_null .OR. ASSOCIATED(norm_drhoa))
352 END IF
353 IF (PRESENT(norm_drhob)) THEN
354 norm_drhob => rho_set%norm_drhob
355 cpassert(my_can_return_null .OR. ASSOCIATED(norm_drhob))
356 END IF
357 IF (PRESENT(rho_1_3)) THEN
358 rho_1_3 => rho_set%rho_1_3
359 cpassert(my_can_return_null .OR. ASSOCIATED(rho_1_3))
360 END IF
361 IF (PRESENT(rhoa_1_3)) THEN
362 rhoa_1_3 => rho_set%rhoa_1_3
363 cpassert(my_can_return_null .OR. ASSOCIATED(rhoa_1_3))
364 END IF
365 IF (PRESENT(rhob_1_3)) THEN
366 rhob_1_3 => rho_set%rhob_1_3
367 cpassert(my_can_return_null .OR. ASSOCIATED(rhob_1_3))
368 END IF
369 IF (PRESENT(tau)) THEN
370 tau => rho_set%tau
371 cpassert(my_can_return_null .OR. ASSOCIATED(tau))
372 END IF
373 IF (PRESENT(tau_a)) THEN
374 tau_a => rho_set%tau_a
375 cpassert(my_can_return_null .OR. ASSOCIATED(tau_a))
376 END IF
377 IF (PRESENT(tau_b)) THEN
378 tau_b => rho_set%tau_b
379 cpassert(my_can_return_null .OR. ASSOCIATED(tau_b))
380 END IF
381 IF (PRESENT(rho_cutoff)) rho_cutoff = rho_set%rho_cutoff
382 IF (PRESENT(drho_cutoff)) drho_cutoff = rho_set%drho_cutoff
383 IF (PRESENT(tau_cutoff)) tau_cutoff = rho_set%tau_cutoff
384 IF (PRESENT(local_bounds)) local_bounds = rho_set%local_bounds
385 END SUBROUTINE xc_rho_set_get
386
387
388! **************************************************************************************************
389!> \brief Shifts association of the requested array to a pw grid
390!> Requires that the corresponding component of rho_set is associated
391!> If owns_data returns TRUE, the caller has to allocate the data later
392!> It is allowed to task for only one component per call.
393!> In case of drho, the array is allocated if not internally available and calculated from drhoa and drhob.
394!> \param rho_set the object you want info about
395!> \param pw_grid ...
396!> \param pw_pool ...
397!> \param owns_data ...
398!> \param rho ...
399!> \param drho ...
400!> \param norm_drho ...
401!> \param rhoa ...
402!> \param rhob ...
403!> \param norm_drhoa ...
404!> \param norm_drhob ...
405!> \param rho_1_3 ...
406!> \param rhoa_1_3 ...
407!> \param rhob_1_3 ...
408!> \param laplace_rho ...
409!> \param laplace_rhoa ...
410!> \param laplace_rhob ...
411!> \param drhoa ...
412!> \param drhob ...
413!> \param tau ...
414!> \param tau_a ...
415!> \param tau_b ...
416! **************************************************************************************************
417 SUBROUTINE xc_rho_set_recover_pw(rho_set, pw_grid, pw_pool, owns_data, rho, drho, norm_drho, &
418 rhoa, rhob, norm_drhoa, norm_drhob, rho_1_3, rhoa_1_3, &
419 rhob_1_3, laplace_rho, laplace_rhoa, laplace_rhob, drhoa, drhob, &
420 tau, tau_a, tau_b)
421 TYPE(xc_rho_set_type) :: rho_set
422 TYPE(pw_r3d_rs_type), DIMENSION(3), OPTIONAL, INTENT(OUT) :: drho, drhoa, drhob
423 TYPE(pw_r3d_rs_type), OPTIONAL, POINTER :: rho, norm_drho, rhoa, rhob, norm_drhoa, &
424 norm_drhob, rho_1_3, rhoa_1_3, &
425 rhob_1_3, laplace_rho, laplace_rhoa, &
426 laplace_rhob, tau, tau_a, tau_b
427 TYPE(pw_grid_type), POINTER, INTENT(IN) :: pw_grid
428 TYPE(pw_pool_type), POINTER, INTENT(IN) :: pw_pool
429 LOGICAL, INTENT(OUT) :: owns_data
430
431 INTEGER :: i
432
433 IF (PRESENT(rho)) THEN
434 NULLIFY (rho)
435 ALLOCATE (rho)
436 CALL xc_rho_set_recover_pw_low(rho, rho_set%rho, pw_grid, pw_pool)
437 NULLIFY (rho_set%rho)
438 owns_data = rho_set%owns%rho
439 END IF
440 IF (PRESENT(drho)) THEN
441 DO i = 1, 3
442 CALL xc_rho_set_recover_pw_low(drho(i), rho_set%drho(i)%array, pw_grid, pw_pool, rho_set%drhoa(i)%array,&
443 & rho_set%drhob(i)%array)
444 END DO
445 owns_data = .true.
446 END IF
447 IF (PRESENT(norm_drho)) THEN
448 NULLIFY (norm_drho)
449 ALLOCATE (norm_drho)
450 CALL xc_rho_set_recover_pw_low(norm_drho, rho_set%norm_drho, pw_grid, pw_pool)
451 NULLIFY (rho_set%norm_drho)
452 owns_data = rho_set%owns%norm_drho
453 END IF
454 IF (PRESENT(rhoa)) THEN
455 NULLIFY (rhoa)
456 ALLOCATE (rhoa)
457 CALL xc_rho_set_recover_pw_low(rhoa, rho_set%rhoa, pw_grid, pw_pool)
458 NULLIFY (rho_set%rhoa)
459 owns_data = rho_set%owns%rho_spin
460 END IF
461 IF (PRESENT(rhob)) THEN
462 NULLIFY (rhob)
463 ALLOCATE (rhob)
464 CALL xc_rho_set_recover_pw_low(rhob, rho_set%rhob, pw_grid, pw_pool)
465 NULLIFY (rho_set%rhob)
466 owns_data = rho_set%owns%rho_spin
467 END IF
468 IF (PRESENT(norm_drhoa)) THEN
469 NULLIFY (norm_drhoa)
470 ALLOCATE (norm_drhoa)
471 CALL xc_rho_set_recover_pw_low(norm_drhoa, rho_set%norm_drhoa, pw_grid, pw_pool)
472 NULLIFY (rho_set%norm_drhoa)
473 owns_data = rho_set%owns%norm_drho_spin
474 END IF
475 IF (PRESENT(norm_drhob)) THEN
476 NULLIFY (norm_drhob)
477 ALLOCATE (norm_drhob)
478 CALL xc_rho_set_recover_pw_low(norm_drhob, rho_set%norm_drhob, pw_grid, pw_pool)
479 NULLIFY (rho_set%norm_drhob)
480 owns_data = rho_set%owns%norm_drho_spin
481 END IF
482 IF (PRESENT(rho_1_3)) THEN
483 NULLIFY (rho_1_3)
484 ALLOCATE (rho_1_3)
485 CALL xc_rho_set_recover_pw_low(rho_1_3, rho_set%rho_1_3, pw_grid, pw_pool)
486 NULLIFY (rho_set%rho_1_3)
487 owns_data = rho_set%owns%rho_1_3
488 END IF
489 IF (PRESENT(rhoa_1_3)) THEN
490 NULLIFY (rhoa_1_3)
491 ALLOCATE (rhoa_1_3)
492 CALL xc_rho_set_recover_pw_low(rhoa_1_3, rho_set%rhoa_1_3, pw_grid, pw_pool)
493 NULLIFY (rho_set%rhoa_1_3)
494 owns_data = rho_set%owns%rho_spin_1_3
495 END IF
496 IF (PRESENT(rhob_1_3)) THEN
497 NULLIFY (rhob_1_3)
498 ALLOCATE (rhob_1_3)
499 CALL xc_rho_set_recover_pw_low(rhob_1_3, rho_set%rhob_1_3, pw_grid, pw_pool)
500 NULLIFY (rho_set%rhob_1_3)
501 owns_data = rho_set%owns%rho_spin_1_3
502 END IF
503 IF (PRESENT(laplace_rho)) THEN
504 NULLIFY (laplace_rho)
505 ALLOCATE (laplace_rho)
506 CALL xc_rho_set_recover_pw_low(laplace_rho, rho_set%laplace_rho, pw_grid, pw_pool)
507 NULLIFY (rho_set%laplace_rho)
508 owns_data = rho_set%owns%laplace_rho
509 END IF
510 IF (PRESENT(laplace_rhoa)) THEN
511 NULLIFY (laplace_rhoa)
512 ALLOCATE (laplace_rhoa)
513 CALL xc_rho_set_recover_pw_low(laplace_rhoa, rho_set%laplace_rhoa, pw_grid, pw_pool)
514 NULLIFY (rho_set%laplace_rhoa)
515 owns_data = rho_set%owns%laplace_rho_spin
516 END IF
517 IF (PRESENT(laplace_rhob)) THEN
518 NULLIFY (laplace_rhob)
519 ALLOCATE (laplace_rhob)
520 CALL xc_rho_set_recover_pw_low(laplace_rhob, rho_set%laplace_rhob, pw_grid, pw_pool)
521 NULLIFY (rho_set%laplace_rhob)
522 owns_data = rho_set%owns%laplace_rho_spin
523 END IF
524 IF (PRESENT(drhoa)) THEN
525 DO i = 1, 3
526 CALL xc_rho_set_recover_pw_low(drhoa(i), rho_set%drhoa(i)%array, pw_grid, pw_pool)
527 END DO
528 owns_data = rho_set%owns%drho_spin
529 END IF
530 IF (PRESENT(drhob)) THEN
531 DO i = 1, 3
532 CALL xc_rho_set_recover_pw_low(drhob(i), rho_set%drhob(i)%array, pw_grid, pw_pool)
533 END DO
534 owns_data = rho_set%owns%drho_spin
535 END IF
536 IF (PRESENT(tau)) THEN
537 NULLIFY (tau)
538 ALLOCATE (tau)
539 CALL xc_rho_set_recover_pw_low(tau, rho_set%tau, pw_grid, pw_pool)
540 NULLIFY (rho_set%tau)
541 owns_data = rho_set%owns%tau
542 END IF
543 IF (PRESENT(tau_a)) THEN
544 NULLIFY (tau_a)
545 ALLOCATE (tau_a)
546 CALL xc_rho_set_recover_pw_low(tau_a, rho_set%tau_a, pw_grid, pw_pool)
547 NULLIFY (rho_set%tau_a)
548 owns_data = rho_set%owns%tau_spin
549 END IF
550 IF (PRESENT(tau_b)) THEN
551 NULLIFY (tau_b)
552 ALLOCATE (tau_b)
553 CALL xc_rho_set_recover_pw_low(tau_b, rho_set%tau_b, pw_grid, pw_pool)
554 NULLIFY (rho_set%tau_b)
555 owns_data = rho_set%owns%tau_spin
556 END IF
557
558 END SUBROUTINE xc_rho_set_recover_pw
559
560 SUBROUTINE xc_rho_set_recover_pw_low(rho_pw, rho, pw_grid, pw_pool, rhoa, rhob)
561 TYPE(pw_r3d_rs_type), INTENT(OUT) :: rho_pw
562 REAL(kind=dp), DIMENSION(:, :, :), POINTER, CONTIGUOUS :: rho
563 TYPE(pw_grid_type), POINTER, INTENT(IN) :: pw_grid
564 TYPE(pw_pool_type), POINTER, INTENT(IN) :: pw_pool
565 REAL(kind=dp), DIMENSION(:, :, :), POINTER, OPTIONAL :: rhoa, rhob
566
567 IF (ASSOCIATED(rho)) THEN
568 CALL rho_pw%create(pw_grid=pw_grid, array_ptr=rho)
569 NULLIFY (rho)
570 ELSE IF (PRESENT(rhoa) .AND. PRESENT(rhob)) THEN
571 CALL pw_pool%create_pw(rho_pw)
572!$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(rho_pw,rhoa,rhob)
573 rho_pw%array(:, :, :) = rhoa(:, :, :) + rhob(:, :, :)
574!$OMP END PARALLEL WORKSHARE
575 ELSE
576 CALL cp_abort(__location__, "Either component or its spin parts (if applicable) "// &
577 "have to be associated in rho_set!")
578 END IF
579
580 END SUBROUTINE xc_rho_set_recover_pw_low
581
582! **************************************************************************************************
583!> \brief cleans (releases) most of the data stored in the rho_set giving back
584!> what it can to the pw_pool
585!> \param rho_set the rho_set to be cleaned
586!> \param pw_pool place to give back 3d arrays,...
587!> \author Fawzi Mohamed
588! **************************************************************************************************
589 SUBROUTINE xc_rho_set_clean(rho_set, pw_pool)
590 TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
591 TYPE(pw_pool_type), POINTER :: pw_pool
592
593 INTEGER :: idir
594
595 IF (rho_set%owns%rho) THEN
596 CALL pw_pool%give_back_cr3d(rho_set%rho)
597 ELSE
598 NULLIFY (rho_set%rho)
599 END IF
600 IF (rho_set%owns%rho_1_3) THEN
601 CALL pw_pool%give_back_cr3d(rho_set%rho_1_3)
602 ELSE
603 NULLIFY (rho_set%rho_1_3)
604 END IF
605 IF (rho_set%owns%drho) THEN
606 DO idir = 1, 3
607 CALL pw_pool%give_back_cr3d(rho_set%drho(idir)%array)
608 END DO
609 ELSE
610 DO idir = 1, 3
611 NULLIFY (rho_set%drho(idir)%array)
612 END DO
613 END IF
614 IF (rho_set%owns%norm_drho) THEN
615 CALL pw_pool%give_back_cr3d(rho_set%norm_drho)
616 ELSE
617 NULLIFY (rho_set%norm_drho)
618 END IF
619 IF (rho_set%owns%laplace_rho) THEN
620 CALL pw_pool%give_back_cr3d(rho_set%laplace_rho)
621 ELSE
622 NULLIFY (rho_set%laplace_rho)
623 END IF
624 IF (rho_set%owns%tau) THEN
625 CALL pw_pool%give_back_cr3d(rho_set%tau)
626 ELSE
627 NULLIFY (rho_set%tau)
628 END IF
629 IF (rho_set%owns%rho_spin) THEN
630 CALL pw_pool%give_back_cr3d(rho_set%rhoa)
631 CALL pw_pool%give_back_cr3d(rho_set%rhob)
632 ELSE
633 NULLIFY (rho_set%rhoa, rho_set%rhob)
634 END IF
635 IF (rho_set%owns%rho_spin_1_3) THEN
636 CALL pw_pool%give_back_cr3d(rho_set%rhoa_1_3)
637 CALL pw_pool%give_back_cr3d(rho_set%rhob_1_3)
638 ELSE
639 NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
640 END IF
641 IF (rho_set%owns%drho_spin) THEN
642 DO idir = 1, 3
643 CALL pw_pool%give_back_cr3d(rho_set%drhoa(idir)%array)
644 CALL pw_pool%give_back_cr3d(rho_set%drhob(idir)%array)
645 END DO
646 ELSE
647 DO idir = 1, 3
648 NULLIFY (rho_set%drhoa(idir)%array, rho_set%drhob(idir)%array)
649 END DO
650 END IF
651 IF (rho_set%owns%laplace_rho_spin) THEN
652 CALL pw_pool%give_back_cr3d(rho_set%laplace_rhoa)
653 CALL pw_pool%give_back_cr3d(rho_set%laplace_rhob)
654 ELSE
655 NULLIFY (rho_set%laplace_rhoa, rho_set%laplace_rhob)
656 END IF
657 IF (rho_set%owns%norm_drho_spin) THEN
658 CALL pw_pool%give_back_cr3d(rho_set%norm_drhoa)
659 CALL pw_pool%give_back_cr3d(rho_set%norm_drhob)
660 ELSE
661 NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
662 END IF
663 IF (rho_set%owns%tau_spin) THEN
664 CALL pw_pool%give_back_cr3d(rho_set%tau_a)
665 CALL pw_pool%give_back_cr3d(rho_set%tau_b)
666 ELSE
667 NULLIFY (rho_set%tau_a, rho_set%tau_b)
668 END IF
669
670 CALL xc_rho_cflags_setall(rho_set%has, .false.)
671 CALL xc_rho_cflags_setall(rho_set%owns, .false.)
672
673 END SUBROUTINE xc_rho_set_clean
674
675! **************************************************************************************************
676!> \brief updates the given rho set with the density given by
677!> rho_r (and rho_g). The rho set will contain the components specified
678!> in needs
679!> \param rho_set the rho_set to update
680!> \param rho_r the new density (in r space)
681!> \param rho_g the new density (in g space, needed for some
682!> derivatives)
683!> \param tau ...
684!> \param needs the components of rho that are needed
685!> \param xc_deriv_method_id ...
686!> \param xc_rho_smooth_id ...
687!> \param pw_pool pool for the allocation of pw and array
688! **************************************************************************************************
689 SUBROUTINE xc_rho_set_update(rho_set, rho_r, rho_g, tau, needs, &
690 xc_deriv_method_id, xc_rho_smooth_id, pw_pool, spinflip)
691 TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
692 TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: rho_r
693 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
694 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER, INTENT(IN) :: tau
695 TYPE(xc_rho_cflags_type), INTENT(in) :: needs
696 INTEGER, INTENT(IN) :: xc_deriv_method_id, xc_rho_smooth_id
697 TYPE(pw_pool_type), POINTER :: pw_pool
698 LOGICAL, OPTIONAL :: spinflip
699
700 REAL(kind=dp), PARAMETER :: f13 = (1.0_dp/3.0_dp)
701
702 INTEGER :: i, idir, ispin, j, k, nspins
703 LOGICAL :: gradient_f, my_rho_g_local, &
704 needs_laplace, needs_rho_g, do_sf
705 REAL(kind=dp) :: rho_cutoff
706 TYPE(pw_r3d_rs_type), DIMENSION(2) :: laplace_rho_r
707 TYPE(pw_r3d_rs_type), DIMENSION(3, 2) :: drho_r
708 TYPE(pw_c1d_gs_type) :: my_rho_g, tmp_g
709 TYPE(pw_r3d_rs_type), DIMENSION(2) :: my_rho_r
710
711 do_sf = .false.
712 IF (PRESENT(spinflip)) do_sf = spinflip
713
714 IF (any(rho_set%local_bounds /= pw_pool%pw_grid%bounds_local)) THEN
715 cpabort("pw_pool cr3d have different size than expected")
716 END IF
717 nspins = SIZE(rho_r)
718 rho_set%local_bounds = rho_r(1)%pw_grid%bounds_local
719 rho_cutoff = 0.5*rho_set%rho_cutoff
720
721 my_rho_g_local = .false.
722 ! some checks
723 SELECT CASE (nspins)
724 CASE (1)
725 IF (.NOT. do_sf) THEN
726 cpassert(.NOT. needs%rho_spin)
727 cpassert(.NOT. needs%drho_spin)
728 cpassert(.NOT. needs%norm_drho_spin)
729 cpassert(.NOT. needs%rho_spin_1_3)
730 cpassert(.NOT. needs%tau_spin)
731 cpassert(.NOT. needs%laplace_rho_spin)
732 ELSE
733 cpassert(.NOT. needs%rho)
734 cpassert(.NOT. needs%drho)
735 cpassert(.NOT. needs%rho_1_3)
736 cpassert(.NOT. needs%tau)
737 cpassert(.NOT. needs%laplace_rho)
738 END IF
739 CASE (2)
740 cpassert(.NOT. needs%rho)
741 cpassert(.NOT. needs%drho)
742 cpassert(.NOT. needs%rho_1_3)
743 cpassert(.NOT. needs%tau)
744 cpassert(.NOT. needs%laplace_rho)
745 CASE default
746 cpabort("Unknown number of spin states")
747 END SELECT
748
749 CALL xc_rho_set_clean(rho_set, pw_pool=pw_pool)
750
751 needs_laplace = (needs%laplace_rho .OR. needs%laplace_rho_spin)
752 gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
753 needs%drho .OR. needs%norm_drho .OR. &
754 needs_laplace)
755 needs_rho_g = ((xc_deriv_method_id == xc_deriv_spline3 .OR. &
756 xc_deriv_method_id == xc_deriv_spline2 .OR. &
757 xc_deriv_method_id == xc_deriv_pw)) .AND. (gradient_f .OR. needs_laplace)
758 IF ((gradient_f .AND. needs_laplace) .AND. &
759 (xc_deriv_method_id /= xc_deriv_pw)) THEN
760 CALL cp_abort(__location__, &
761 "MGGA functionals that require the Laplacian are "// &
762 "only compatible with 'XC_DERIV PW' and 'XC_SMOOTH_RHO NONE'")
763 END IF
764
765 IF (needs_rho_g) THEN
766 CALL pw_pool%create_pw(tmp_g)
767 END IF
768 DO ispin = 1, nspins
769 CALL pw_pool%create_pw(my_rho_r(ispin))
770 ! introduce a smoothing kernel on the density
771 IF (xc_rho_smooth_id == xc_rho_no_smooth) THEN
772 IF (needs_rho_g) THEN
773 IF (ASSOCIATED(rho_g)) THEN
774 my_rho_g_local = .false.
775 my_rho_g = rho_g(ispin)
776 END IF
777 END IF
778
779 CALL pw_copy(rho_r(ispin), my_rho_r(ispin))
780 ELSE
781 CALL xc_pw_smooth(rho_r(ispin), my_rho_r(ispin), xc_rho_smooth_id)
782 END IF
783
784 IF (gradient_f) THEN ! calculate the grad of rho
785 ! normally when you need the gradient you need the whole gradient
786 ! (for the partial integration)
787 ! deriv rho
788 DO idir = 1, 3
789 CALL pw_pool%create_pw(drho_r(idir, ispin))
790 END DO
791 IF (needs_rho_g) THEN
792 IF (.NOT. ASSOCIATED(my_rho_g%pw_grid)) THEN
793 my_rho_g_local = .true.
794 CALL pw_pool%create_pw(my_rho_g)
795 CALL pw_transfer(my_rho_r(ispin), my_rho_g)
796 END IF
797 IF (.NOT. my_rho_g_local .AND. (xc_deriv_method_id == xc_deriv_spline2 .OR. &
798 xc_deriv_method_id == xc_deriv_spline3)) THEN
799 CALL pw_pool%create_pw(my_rho_g)
800 my_rho_g_local = .true.
801 CALL pw_copy(rho_g(ispin), my_rho_g)
802 END IF
803 END IF
804 IF (needs%laplace_rho .OR. needs%laplace_rho_spin) THEN
805 CALL pw_pool%create_pw(laplace_rho_r(ispin))
806 CALL xc_pw_laplace(my_rho_g, pw_pool, xc_deriv_method_id, laplace_rho_r(ispin), tmp_g=tmp_g)
807 END IF
808 CALL xc_pw_gradient(my_rho_r(ispin), my_rho_g, tmp_g, drho_r(:, ispin), xc_deriv_method_id)
809
810 IF (needs_rho_g) THEN
811 IF (my_rho_g_local) THEN
812 my_rho_g_local = .false.
813 CALL pw_pool%give_back_pw(my_rho_g)
814 END IF
815 END IF
816
817 IF (xc_deriv_method_id /= xc_deriv_pw) THEN
818 CALL pw_spline_scale_deriv(drho_r(:, ispin))
819 END IF
820
821 END IF
822
823 END DO
824
825 IF (ASSOCIATED(tmp_g%pw_grid)) THEN
826 CALL pw_pool%give_back_pw(tmp_g)
827 END IF
828
829 SELECT CASE (nspins)
830 CASE (1)
831 IF (.NOT. do_sf) THEN
832 IF (needs%rho_1_3) THEN
833 CALL pw_pool%create_cr3d(rho_set%rho_1_3)
834!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,my_rho_r)
835 DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
836 DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
837 DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
838 rho_set%rho_1_3(i, j, k) = max(my_rho_r(1)%array(i, j, k), 0.0_dp)**f13
839 END DO
840 END DO
841 END DO
842 rho_set%owns%rho_1_3 = .true.
843 rho_set%has%rho_1_3 = .true.
844 END IF
845 IF (needs%rho) THEN
846 rho_set%rho => my_rho_r(1)%array
847 NULLIFY (my_rho_r(1)%array)
848 rho_set%owns%rho = .true.
849 rho_set%has%rho = .true.
850 END IF
851 IF (needs%norm_drho) THEN
852 CALL pw_pool%create_cr3d(rho_set%norm_drho)
853!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
854 DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
855 DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
856 DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
857 rho_set%norm_drho(i, j, k) = sqrt( &
858 drho_r(1, 1)%array(i, j, k)**2 + &
859 drho_r(2, 1)%array(i, j, k)**2 + &
860 drho_r(3, 1)%array(i, j, k)**2)
861 END DO
862 END DO
863 END DO
864 rho_set%owns%norm_drho = .true.
865 rho_set%has%norm_drho = .true.
866 END IF
867 IF (needs%laplace_rho) THEN
868 rho_set%laplace_rho => laplace_rho_r(1)%array
869 NULLIFY (laplace_rho_r(1)%array)
870 rho_set%owns%laplace_rho = .true.
871 rho_set%has%laplace_rho = .true.
872 END IF
873
874 IF (needs%drho) THEN
875 DO idir = 1, 3
876 rho_set%drho(idir)%array => drho_r(idir, 1)%array
877 NULLIFY (drho_r(idir, 1)%array)
878 END DO
879 rho_set%owns%drho = .true.
880 rho_set%has%drho = .true.
881 END IF
882 ELSE
883 IF (needs%norm_drho) THEN
884 CALL pw_pool%create_cr3d(rho_set%norm_drho)
885!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
886 DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
887 DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
888 DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
889 rho_set%norm_drho(i, j, k) = sqrt( &
890 drho_r(1, 1)%array(i, j, k)**2 + &
891 drho_r(2, 1)%array(i, j, k)**2 + &
892 drho_r(3, 1)%array(i, j, k)**2)
893 END DO
894 END DO
895 END DO
896 rho_set%owns%norm_drho = .true.
897 rho_set%has%norm_drho = .true.
898 END IF
899 IF (needs%rho_spin) THEN
900
901 rho_set%rhoa => my_rho_r(1)%array
902 NULLIFY (my_rho_r(1)%array)
903
904 rho_set%owns%rho_spin = .true.
905 rho_set%has%rho_spin = .true.
906 END IF
907 IF (needs%norm_drho_spin) THEN
908 CALL pw_pool%create_cr3d(rho_set%norm_drhoa)
909!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
910 DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
911 DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
912 DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
913 rho_set%norm_drhoa(i, j, k) = sqrt( &
914 drho_r(1, 1)%array(i, j, k)**2 + &
915 drho_r(2, 1)%array(i, j, k)**2 + &
916 drho_r(3, 1)%array(i, j, k)**2)
917 END DO
918 END DO
919 END DO
920 rho_set%owns%norm_drho_spin = .true.
921 rho_set%has%norm_drho_spin = .true.
922 END IF
923 IF (needs%laplace_rho_spin) THEN
924 rho_set%laplace_rhoa => laplace_rho_r(1)%array
925 NULLIFY (laplace_rho_r(1)%array)
926
927 rho_set%owns%laplace_rho_spin = .true.
928 rho_set%has%laplace_rho_spin = .true.
929 END IF
930 IF (needs%drho_spin) THEN
931 DO idir = 1, 3
932 rho_set%drhoa(idir)%array => drho_r(idir, 1)%array
933 NULLIFY (drho_r(idir, 1)%array)
934 END DO
935 rho_set%owns%drho_spin = .true.
936 rho_set%has%drho_spin = .true.
937 END IF
938 END IF
939 CASE (2)
940 IF (needs%rho_spin_1_3) THEN
941 CALL pw_pool%create_cr3d(rho_set%rhoa_1_3)
942 !assume that the bounds are the same?
943!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,my_rho_r)
944 DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
945 DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
946 DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
947 rho_set%rhoa_1_3(i, j, k) = max(my_rho_r(1)%array(i, j, k), 0.0_dp)**f13
948 END DO
949 END DO
950 END DO
951 CALL pw_pool%create_cr3d(rho_set%rhob_1_3)
952 !assume that the bounds are the same?
953!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,my_rho_r)
954 DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
955 DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
956 DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
957 rho_set%rhob_1_3(i, j, k) = max(my_rho_r(2)%array(i, j, k), 0.0_dp)**f13
958 END DO
959 END DO
960 END DO
961 rho_set%owns%rho_spin_1_3 = .true.
962 rho_set%has%rho_spin_1_3 = .true.
963 END IF
964 IF (needs%norm_drho) THEN
965
966 CALL pw_pool%create_cr3d(rho_set%norm_drho)
967!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
968 DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
969 DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
970 DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
971 rho_set%norm_drho(i, j, k) = sqrt( &
972 (drho_r(1, 1)%array(i, j, k) + drho_r(1, 2)%array(i, j, k))**2 + &
973 (drho_r(2, 1)%array(i, j, k) + drho_r(2, 2)%array(i, j, k))**2 + &
974 (drho_r(3, 1)%array(i, j, k) + drho_r(3, 2)%array(i, j, k))**2)
975 END DO
976 END DO
977 END DO
978
979 rho_set%owns%norm_drho = .true.
980 rho_set%has%norm_drho = .true.
981 END IF
982 IF (needs%rho_spin) THEN
983
984 rho_set%rhoa => my_rho_r(1)%array
985 NULLIFY (my_rho_r(1)%array)
986
987 rho_set%rhob => my_rho_r(2)%array
988 NULLIFY (my_rho_r(2)%array)
989
990 rho_set%owns%rho_spin = .true.
991 rho_set%has%rho_spin = .true.
992 END IF
993 IF (needs%norm_drho_spin) THEN
994
995 CALL pw_pool%create_cr3d(rho_set%norm_drhoa)
996!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
997 DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
998 DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
999 DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
1000 rho_set%norm_drhoa(i, j, k) = sqrt( &
1001 drho_r(1, 1)%array(i, j, k)**2 + &
1002 drho_r(2, 1)%array(i, j, k)**2 + &
1003 drho_r(3, 1)%array(i, j, k)**2)
1004 END DO
1005 END DO
1006 END DO
1007
1008 CALL pw_pool%create_cr3d(rho_set%norm_drhob)
1009 rho_set%owns%norm_drho_spin = .true.
1010!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
1011 DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
1012 DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
1013 DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
1014 rho_set%norm_drhob(i, j, k) = sqrt( &
1015 drho_r(1, 2)%array(i, j, k)**2 + &
1016 drho_r(2, 2)%array(i, j, k)**2 + &
1017 drho_r(3, 2)%array(i, j, k)**2)
1018 END DO
1019 END DO
1020 END DO
1021
1022 rho_set%owns%norm_drho_spin = .true.
1023 rho_set%has%norm_drho_spin = .true.
1024 END IF
1025 IF (needs%laplace_rho_spin) THEN
1026 rho_set%laplace_rhoa => laplace_rho_r(1)%array
1027 NULLIFY (laplace_rho_r(1)%array)
1028
1029 rho_set%laplace_rhob => laplace_rho_r(2)%array
1030 NULLIFY (laplace_rho_r(2)%array)
1031
1032 rho_set%owns%laplace_rho_spin = .true.
1033 rho_set%has%laplace_rho_spin = .true.
1034 END IF
1035 IF (needs%drho_spin) THEN
1036 DO idir = 1, 3
1037 rho_set%drhoa(idir)%array => drho_r(idir, 1)%array
1038 NULLIFY (drho_r(idir, 1)%array)
1039 rho_set%drhob(idir)%array => drho_r(idir, 2)%array
1040 NULLIFY (drho_r(idir, 2)%array)
1041 END DO
1042 rho_set%owns%drho_spin = .true.
1043 rho_set%has%drho_spin = .true.
1044 END IF
1045 END SELECT
1046 ! post cleanup
1047 DO ispin = 1, nspins
1048 IF (needs%laplace_rho .OR. needs%laplace_rho_spin) THEN
1049 CALL pw_pool%give_back_pw(laplace_rho_r(ispin))
1050 END IF
1051 DO idir = 1, 3
1052 CALL pw_pool%give_back_pw(drho_r(idir, ispin))
1053 END DO
1054 END DO
1055 DO ispin = 1, nspins
1056 CALL pw_pool%give_back_pw(my_rho_r(ispin))
1057 END DO
1058
1059 ! tau part
1060 IF (needs%tau .OR. needs%tau_spin) THEN
1061 cpassert(ASSOCIATED(tau))
1062 DO ispin = 1, nspins
1063 cpassert(ASSOCIATED(tau(ispin)%array))
1064 END DO
1065 END IF
1066 IF (needs%tau) THEN
1067 rho_set%tau => tau(1)%array
1068 rho_set%owns%tau = .false.
1069 rho_set%has%tau = .true.
1070 END IF
1071 IF (needs%tau_spin) THEN
1072 rho_set%tau_a => tau(1)%array
1073 rho_set%tau_b => tau(2)%array
1074 rho_set%owns%tau_spin = .false.
1075 rho_set%has%tau_spin = .true.
1076 END IF
1077
1078 cpassert(xc_rho_cflags_equal(rho_set%has, needs))
1079
1080 END SUBROUTINE xc_rho_set_update
1081
1082 END MODULE xc_rho_set_types
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_transfer
Definition pw_methods.F:304
xc_util::xc_pw_laplace
Definition xc_util.F:46
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition cp_array_utils.F:18
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
pw_grid_types
Definition pw_grid_types.F:13
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_spline_utils
different utils that are useful to manipulate splines on the regular grid of a pw
Definition pw_spline_utils.F:18
pw_spline_utils::pw_spline_scale_deriv
subroutine, public pw_spline_scale_deriv(deriv_vals_r, transpose, scale)
rescales the derivatives from gridspacing=1 to the real derivatives
Definition pw_spline_utils.F:275
pw_types
Definition pw_types.F:24
xc_input_constants
input constants for xc
Definition xc_input_constants.F:10
xc_input_constants::xc_deriv_pw
integer, parameter, public xc_deriv_pw
Definition xc_input_constants.F:14
xc_input_constants::xc_deriv_spline2
integer, parameter, public xc_deriv_spline2
Definition xc_input_constants.F:14
xc_input_constants::xc_deriv_spline3
integer, parameter, public xc_deriv_spline3
Definition xc_input_constants.F:14
xc_input_constants::xc_rho_no_smooth
integer, parameter, public xc_rho_no_smooth
Definition xc_input_constants.F:23
xc_rho_cflags_types
contains the structure
Definition xc_rho_cflags_types.F:14
xc_rho_cflags_types::xc_rho_cflags_equal
elemental logical function, public xc_rho_cflags_equal(cflags1, cflags2)
return true if the two cflags are equal
Definition xc_rho_cflags_types.F:87
xc_rho_cflags_types::xc_rho_cflags_setall
elemental subroutine, public xc_rho_cflags_setall(cflags, value)
sets all the flags to the given value
Definition xc_rho_cflags_types.F:63
xc_rho_set_types
contains the structure
Definition xc_rho_set_types.F:14
xc_rho_set_types::xc_rho_set_create
subroutine, public xc_rho_set_create(rho_set, local_bounds, rho_cutoff, drho_cutoff, tau_cutoff)
allocates and does (minimal) initialization of a rho_set
Definition xc_rho_set_types.F:111
xc_rho_set_types::xc_rho_set_release
subroutine, public xc_rho_set_release(rho_set, pw_pool)
releases the given rho_set
Definition xc_rho_set_types.F:129
xc_rho_set_types::xc_rho_set_recover_pw
subroutine, public xc_rho_set_recover_pw(rho_set, pw_grid, pw_pool, owns_data, 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)
Shifts association of the requested array to a pw grid Requires that the corresponding component of r...
Definition xc_rho_set_types.F:421
xc_rho_set_types::xc_rho_set_update
subroutine, public xc_rho_set_update(rho_set, rho_r, rho_g, tau, needs, xc_deriv_method_id, xc_rho_smooth_id, pw_pool, spinflip)
updates the given rho set with the density given by rho_r (and rho_g). The rho set will contain the c...
Definition xc_rho_set_types.F:691
xc_rho_set_types::xc_rho_set_get
subroutine, public xc_rho_set_get(rho_set, can_return_null, 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, rho_cutoff, drho_cutoff, tau_cutoff, tau, tau_a, tau_b, local_bounds)
returns the various attributes of rho_set
Definition xc_rho_set_types.F:281
xc_util
contains utility functions for the xc package
Definition xc_util.F:14
xc_util::xc_pw_smooth
subroutine, public xc_pw_smooth(pw_in, pw_out, xc_smooth_id)
...
Definition xc_util.F:73
xc_util::xc_pw_gradient
subroutine, public xc_pw_gradient(pw_r, pw_g, tmp_g, gradient, xc_deriv_method_id)
...
Definition xc_util.F:123
cp_array_utils::cp_3d_r_cp_type
represent a pointer to a contiguous 3d array
Definition cp_array_utils.F:149
pw_grid_types::pw_grid_type
Definition pw_grid_types.F:53
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78