Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 367154 : 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 367154 : IF (PRESENT(rho_cutoff)) rho_set%rho_cutoff = rho_cutoff
116 367154 : IF (PRESENT(drho_cutoff)) rho_set%drho_cutoff = drho_cutoff
117 367154 : IF (PRESENT(tau_cutoff)) rho_set%tau_cutoff = tau_cutoff
118 3671540 : rho_set%local_bounds = local_bounds
119 367154 : CALL xc_rho_cflags_setall(rho_set%owns, .TRUE.)
120 367154 : CALL xc_rho_cflags_setall(rho_set%has, .FALSE.)
121 367154 : 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 367154 : 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 367154 : IF (PRESENT(pw_pool)) THEN
135 157871 : IF (ASSOCIATED(pw_pool)) THEN
136 157871 : 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 1468616 : rho_set%local_bounds(1, :) = -HUGE(0) ! we want to crash...
143 1468616 : rho_set%local_bounds(1, :) = HUGE(0)
144 367154 : IF (rho_set%owns%rho .AND. ASSOCIATED(rho_set%rho)) THEN
145 183934 : DEALLOCATE (rho_set%rho)
146 : ELSE
147 183220 : NULLIFY (rho_set%rho)
148 : END IF
149 367154 : IF (rho_set%owns%rho_spin) THEN
150 182065 : IF (ASSOCIATED(rho_set%rhoa)) THEN
151 25349 : DEALLOCATE (rho_set%rhoa)
152 : END IF
153 182065 : IF (ASSOCIATED(rho_set%rhob)) THEN
154 25241 : DEALLOCATE (rho_set%rhob)
155 : END IF
156 : ELSE
157 185089 : NULLIFY (rho_set%rhoa, rho_set%rhob)
158 : END IF
159 367154 : IF (rho_set%owns%rho_1_3 .AND. ASSOCIATED(rho_set%rho_1_3)) THEN
160 5725 : DEALLOCATE (rho_set%rho_1_3)
161 : ELSE
162 361429 : NULLIFY (rho_set%rho_1_3)
163 : END IF
164 367154 : IF (rho_set%owns%rho_spin) THEN
165 182065 : IF (ASSOCIATED(rho_set%rhoa_1_3)) THEN
166 2562 : DEALLOCATE (rho_set%rhoa_1_3)
167 : END IF
168 182065 : IF (ASSOCIATED(rho_set%rhob_1_3)) THEN
169 2562 : DEALLOCATE (rho_set%rhob_1_3)
170 : END IF
171 : ELSE
172 185089 : NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
173 : END IF
174 367154 : IF (rho_set%owns%drho) THEN
175 767532 : DO i = 1, 3
176 767532 : IF (ASSOCIATED(rho_set%drho(i)%array)) THEN
177 326718 : DEALLOCATE (rho_set%drho(i)%array)
178 : END IF
179 : END DO
180 : ELSE
181 701084 : DO i = 1, 3
182 701084 : NULLIFY (rho_set%drho(i)%array)
183 : END DO
184 : END IF
185 367154 : IF (rho_set%owns%drho_spin) THEN
186 716364 : DO i = 1, 3
187 537273 : IF (ASSOCIATED(rho_set%drhoa(i)%array)) THEN
188 46434 : DEALLOCATE (rho_set%drhoa(i)%array)
189 : END IF
190 716364 : IF (ASSOCIATED(rho_set%drhob(i)%array)) THEN
191 46266 : DEALLOCATE (rho_set%drhob(i)%array)
192 : END IF
193 : END DO
194 : ELSE
195 752252 : DO i = 1, 3
196 752252 : NULLIFY (rho_set%drhoa(i)%array, rho_set%drhob(i)%array)
197 : END DO
198 : END IF
199 367154 : IF (rho_set%owns%laplace_rho .AND. ASSOCIATED(rho_set%laplace_rho)) THEN
200 302 : DEALLOCATE (rho_set%laplace_rho)
201 : ELSE
202 366852 : NULLIFY (rho_set%laplace_rho)
203 : END IF
204 :
205 367154 : IF (rho_set%owns%norm_drho .AND. ASSOCIATED(rho_set%norm_drho)) THEN
206 146251 : DEALLOCATE (rho_set%norm_drho)
207 : ELSE
208 220903 : NULLIFY (rho_set%norm_drho)
209 : END IF
210 367154 : IF (rho_set%owns%laplace_rho_spin) THEN
211 175525 : IF (ASSOCIATED(rho_set%laplace_rhoa)) THEN
212 160 : DEALLOCATE (rho_set%laplace_rhoa)
213 : END IF
214 175525 : IF (ASSOCIATED(rho_set%laplace_rhob)) THEN
215 160 : DEALLOCATE (rho_set%laplace_rhob)
216 : END IF
217 : ELSE
218 191629 : NULLIFY (rho_set%laplace_rhoa, rho_set%laplace_rhob)
219 : END IF
220 :
221 367154 : IF (rho_set%owns%norm_drho_spin) THEN
222 179091 : IF (ASSOCIATED(rho_set%norm_drhoa)) THEN
223 16357 : DEALLOCATE (rho_set%norm_drhoa)
224 : END IF
225 179091 : IF (ASSOCIATED(rho_set%norm_drhob)) THEN
226 16301 : DEALLOCATE (rho_set%norm_drhob)
227 : END IF
228 : ELSE
229 188063 : NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
230 : END IF
231 367154 : IF (rho_set%owns%tau .AND. ASSOCIATED(rho_set%tau)) THEN
232 2728 : DEALLOCATE (rho_set%tau)
233 : ELSE
234 364426 : NULLIFY (rho_set%tau)
235 : END IF
236 367154 : IF (rho_set%owns%tau_spin) THEN
237 175365 : IF (ASSOCIATED(rho_set%tau_a)) THEN
238 34 : DEALLOCATE (rho_set%tau_a)
239 : END IF
240 175365 : IF (ASSOCIATED(rho_set%tau_b)) THEN
241 34 : DEALLOCATE (rho_set%tau_b)
242 : END IF
243 : ELSE
244 191789 : NULLIFY (rho_set%tau_a, rho_set%tau_b)
245 : END IF
246 367154 : 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 1412862 : 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 1412862 : my_can_return_null = .FALSE.
301 1412862 : IF (PRESENT(can_return_null)) my_can_return_null = can_return_null
302 :
303 1412862 : IF (PRESENT(rho)) THEN
304 367981 : rho => rho_set%rho
305 367981 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rho))
306 : END IF
307 1412862 : IF (PRESENT(drho)) THEN
308 956188 : DO i = 1, 3
309 717141 : drho(i)%array => rho_set%drho(i)%array
310 956188 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rho_set%drho(i)%array))
311 : END DO
312 : END IF
313 1412862 : IF (PRESENT(norm_drho)) THEN
314 508195 : norm_drho => rho_set%norm_drho
315 508195 : CPASSERT(my_can_return_null .OR. ASSOCIATED(norm_drho))
316 : END IF
317 1412862 : IF (PRESENT(laplace_rho)) THEN
318 18910 : laplace_rho => rho_set%laplace_rho
319 18910 : CPASSERT(my_can_return_null .OR. ASSOCIATED(laplace_rho))
320 : END IF
321 1412862 : IF (PRESENT(rhoa)) THEN
322 190620 : rhoa => rho_set%rhoa
323 190620 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rhoa))
324 : END IF
325 1412862 : IF (PRESENT(rhob)) THEN
326 190008 : rhob => rho_set%rhob
327 190008 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rhob))
328 : END IF
329 1412862 : IF (PRESENT(drhoa)) THEN
330 775636 : DO i = 1, 3
331 581727 : drhoa(i)%array => rho_set%drhoa(i)%array
332 775636 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rho_set%drhoa(i)%array))
333 : END DO
334 : END IF
335 1412862 : IF (PRESENT(drhob)) THEN
336 773796 : DO i = 1, 3
337 580347 : drhob(i)%array => rho_set%drhob(i)%array
338 773796 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rho_set%drhob(i)%array))
339 : END DO
340 : END IF
341 1412862 : IF (PRESENT(laplace_rhoa)) THEN
342 3850 : laplace_rhoa => rho_set%laplace_rhoa
343 3850 : CPASSERT(my_can_return_null .OR. ASSOCIATED(laplace_rhoa))
344 : END IF
345 1412862 : IF (PRESENT(laplace_rhob)) THEN
346 3850 : laplace_rhob => rho_set%laplace_rhob
347 3850 : CPASSERT(my_can_return_null .OR. ASSOCIATED(laplace_rhob))
348 : END IF
349 1412862 : IF (PRESENT(norm_drhoa)) THEN
350 257883 : norm_drhoa => rho_set%norm_drhoa
351 257883 : CPASSERT(my_can_return_null .OR. ASSOCIATED(norm_drhoa))
352 : END IF
353 1412862 : IF (PRESENT(norm_drhob)) THEN
354 257883 : norm_drhob => rho_set%norm_drhob
355 257883 : CPASSERT(my_can_return_null .OR. ASSOCIATED(norm_drhob))
356 : END IF
357 1412862 : IF (PRESENT(rho_1_3)) THEN
358 22695 : rho_1_3 => rho_set%rho_1_3
359 22695 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rho_1_3))
360 : END IF
361 1412862 : IF (PRESENT(rhoa_1_3)) THEN
362 3380 : rhoa_1_3 => rho_set%rhoa_1_3
363 3380 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rhoa_1_3))
364 : END IF
365 1412862 : IF (PRESENT(rhob_1_3)) THEN
366 3380 : rhob_1_3 => rho_set%rhob_1_3
367 3380 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rhob_1_3))
368 : END IF
369 1412862 : IF (PRESENT(tau)) THEN
370 22716 : tau => rho_set%tau
371 22716 : CPASSERT(my_can_return_null .OR. ASSOCIATED(tau))
372 : END IF
373 1412862 : IF (PRESENT(tau_a)) THEN
374 4214 : tau_a => rho_set%tau_a
375 4214 : CPASSERT(my_can_return_null .OR. ASSOCIATED(tau_a))
376 : END IF
377 1412862 : IF (PRESENT(tau_b)) THEN
378 4214 : tau_b => rho_set%tau_b
379 4214 : CPASSERT(my_can_return_null .OR. ASSOCIATED(tau_b))
380 : END IF
381 1412862 : IF (PRESENT(rho_cutoff)) rho_cutoff = rho_set%rho_cutoff
382 1412862 : IF (PRESENT(drho_cutoff)) drho_cutoff = rho_set%drho_cutoff
383 1412862 : IF (PRESENT(tau_cutoff)) tau_cutoff = rho_set%tau_cutoff
384 3628032 : IF (PRESENT(local_bounds)) local_bounds = rho_set%local_bounds
385 1412862 : 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 585765 : 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 585765 : $:recover(variable)
465 : #:endfor
466 :
467 117153 : END SUBROUTINE xc_rho_set_recover_pw
468 :
469 351459 : 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 351459 : IF (ASSOCIATED(rho)) THEN
477 300483 : CALL rho_pw%create(pw_grid=pw_grid, array_ptr=rho)
478 300483 : NULLIFY (rho)
479 50976 : ELSE IF (PRESENT(rhoa) .AND. PRESENT(rhob)) THEN
480 50976 : CALL pw_pool%create_pw(rho_pw)
481 50976 : !$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 351459 : 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 349660 : 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 349660 : IF (rho_set%owns%rho) THEN
505 314569 : CALL pw_pool%give_back_cr3d(rho_set%rho)
506 : ELSE
507 35091 : NULLIFY (rho_set%rho)
508 : END IF
509 349660 : IF (rho_set%owns%rho_1_3) THEN
510 200561 : CALL pw_pool%give_back_cr3d(rho_set%rho_1_3)
511 : ELSE
512 149099 : NULLIFY (rho_set%rho_1_3)
513 : END IF
514 349660 : IF (rho_set%owns%drho) THEN
515 1025896 : DO idir = 1, 3
516 1025896 : CALL pw_pool%give_back_cr3d(rho_set%drho(idir)%array)
517 : END DO
518 : ELSE
519 372744 : DO idir = 1, 3
520 372744 : NULLIFY (rho_set%drho(idir)%array)
521 : END DO
522 : END IF
523 349660 : IF (rho_set%owns%norm_drho) THEN
524 279580 : CALL pw_pool%give_back_cr3d(rho_set%norm_drho)
525 : ELSE
526 70080 : NULLIFY (rho_set%norm_drho)
527 : END IF
528 349660 : IF (rho_set%owns%laplace_rho) THEN
529 192591 : CALL pw_pool%give_back_cr3d(rho_set%laplace_rho)
530 : ELSE
531 157069 : NULLIFY (rho_set%laplace_rho)
532 : END IF
533 349660 : IF (rho_set%owns%tau) THEN
534 191789 : CALL pw_pool%give_back_cr3d(rho_set%tau)
535 : ELSE
536 157871 : NULLIFY (rho_set%tau)
537 : END IF
538 349660 : IF (rho_set%owns%rho_spin) THEN
539 226880 : CALL pw_pool%give_back_cr3d(rho_set%rhoa)
540 226880 : CALL pw_pool%give_back_cr3d(rho_set%rhob)
541 : ELSE
542 122780 : NULLIFY (rho_set%rhoa, rho_set%rhob)
543 : END IF
544 349660 : IF (rho_set%owns%rho_spin_1_3) THEN
545 193451 : CALL pw_pool%give_back_cr3d(rho_set%rhoa_1_3)
546 193451 : CALL pw_pool%give_back_cr3d(rho_set%rhob_1_3)
547 : ELSE
548 156209 : NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
549 : END IF
550 349660 : IF (rho_set%owns%drho_spin) THEN
551 838828 : DO idir = 1, 3
552 629121 : CALL pw_pool%give_back_cr3d(rho_set%drhoa(idir)%array)
553 838828 : CALL pw_pool%give_back_cr3d(rho_set%drhob(idir)%array)
554 : END DO
555 : ELSE
556 559812 : DO idir = 1, 3
557 559812 : NULLIFY (rho_set%drhoa(idir)%array, rho_set%drhob(idir)%array)
558 : END DO
559 : END IF
560 349660 : IF (rho_set%owns%laplace_rho_spin) THEN
561 192181 : CALL pw_pool%give_back_cr3d(rho_set%laplace_rhoa)
562 192181 : CALL pw_pool%give_back_cr3d(rho_set%laplace_rhob)
563 : ELSE
564 157479 : NULLIFY (rho_set%laplace_rhoa, rho_set%laplace_rhob)
565 : END IF
566 349660 : IF (rho_set%owns%norm_drho_spin) THEN
567 212901 : CALL pw_pool%give_back_cr3d(rho_set%norm_drhoa)
568 212901 : CALL pw_pool%give_back_cr3d(rho_set%norm_drhob)
569 : ELSE
570 136759 : NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
571 : END IF
572 349660 : IF (rho_set%owns%tau_spin) THEN
573 191789 : CALL pw_pool%give_back_cr3d(rho_set%tau_a)
574 191789 : CALL pw_pool%give_back_cr3d(rho_set%tau_b)
575 : ELSE
576 157871 : NULLIFY (rho_set%tau_a, rho_set%tau_b)
577 : END IF
578 :
579 349660 : CALL xc_rho_cflags_setall(rho_set%has, .FALSE.)
580 349660 : CALL xc_rho_cflags_setall(rho_set%owns, .FALSE.)
581 :
582 349660 : 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 191789 : 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 575367 : TYPE(pw_r3d_rs_type), DIMENSION(2) :: laplace_rho_r
616 1726101 : TYPE(pw_r3d_rs_type), DIMENSION(3, 2) :: drho_r
617 : TYPE(pw_c1d_gs_type) :: my_rho_g, tmp_g
618 575367 : TYPE(pw_r3d_rs_type), DIMENSION(2) :: my_rho_r
619 :
620 191789 : do_sf = .FALSE.
621 26944 : IF (PRESENT(spinflip)) do_sf = spinflip
622 :
623 1917890 : IF (ANY(rho_set%local_bounds /= pw_pool%pw_grid%bounds_local)) THEN
624 0 : CPABORT("pw_pool cr3d have different size than expected")
625 : END IF
626 191789 : nspins = SIZE(rho_r)
627 1917890 : rho_set%local_bounds = rho_r(1)%pw_grid%bounds_local
628 191789 : rho_cutoff = 0.5_dp*rho_set%rho_cutoff
629 :
630 191789 : my_rho_g_local = .FALSE.
631 : ! some checks
632 150106 : SELECT CASE (nspins)
633 : CASE (1)
634 150106 : IF (.NOT. do_sf) THEN
635 149998 : CPASSERT(.NOT. needs%rho_spin)
636 149998 : CPASSERT(.NOT. needs%drho_spin)
637 149998 : CPASSERT(.NOT. needs%norm_drho_spin)
638 149998 : CPASSERT(.NOT. needs%rho_spin_1_3)
639 149998 : CPASSERT(.NOT. needs%tau_spin)
640 149998 : CPASSERT(.NOT. needs%laplace_rho_spin)
641 : ELSE
642 108 : CPASSERT(.NOT. needs%rho)
643 108 : CPASSERT(.NOT. needs%drho)
644 108 : CPASSERT(.NOT. needs%rho_1_3)
645 108 : CPASSERT(.NOT. needs%tau)
646 108 : CPASSERT(.NOT. needs%laplace_rho)
647 : END IF
648 : CASE (2)
649 41683 : CPASSERT(.NOT. needs%rho)
650 41683 : CPASSERT(.NOT. needs%drho)
651 41683 : CPASSERT(.NOT. needs%rho_1_3)
652 41683 : CPASSERT(.NOT. needs%tau)
653 41683 : CPASSERT(.NOT. needs%laplace_rho)
654 : CASE default
655 191789 : CPABORT("Unknown number of spin states")
656 : END SELECT
657 :
658 191789 : CALL xc_rho_set_clean(rho_set, pw_pool=pw_pool)
659 :
660 191789 : needs_laplace = (needs%laplace_rho .OR. needs%laplace_rho_spin)
661 : gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
662 : needs%drho .OR. needs%norm_drho .OR. &
663 191789 : needs_laplace)
664 : needs_rho_g = ((xc_deriv_method_id == xc_deriv_spline3 .OR. &
665 : xc_deriv_method_id == xc_deriv_spline2 .OR. &
666 191789 : xc_deriv_method_id == xc_deriv_pw)) .AND. (gradient_f .OR. needs_laplace)
667 109937 : IF ((gradient_f .AND. needs_laplace) .AND. &
668 : (xc_deriv_method_id /= xc_deriv_pw)) THEN
669 : CALL cp_abort(__LOCATION__, &
670 : "MGGA functionals that require the Laplacian are "// &
671 0 : "only compatible with 'XC_DERIV PW' and 'XC_SMOOTH_RHO NONE'")
672 : END IF
673 :
674 191789 : IF (needs_rho_g) THEN
675 107963 : CALL pw_pool%create_pw(tmp_g)
676 : END IF
677 425261 : DO ispin = 1, nspins
678 233472 : CALL pw_pool%create_pw(my_rho_r(ispin))
679 : ! introduce a smoothing kernel on the density
680 233472 : IF (xc_rho_smooth_id == xc_rho_no_smooth) THEN
681 232896 : IF (needs_rho_g) THEN
682 132145 : IF (ASSOCIATED(rho_g)) THEN
683 105417 : my_rho_g_local = .FALSE.
684 105417 : my_rho_g = rho_g(ispin)
685 : END IF
686 : END IF
687 :
688 232896 : CALL pw_copy(rho_r(ispin), my_rho_r(ispin))
689 : ELSE
690 576 : CALL xc_pw_smooth(rho_r(ispin), my_rho_r(ispin), xc_rho_smooth_id)
691 : END IF
692 :
693 425261 : IF (gradient_f) THEN ! calculate the grad of rho
694 : ! normally when you need the gradient you need the whole gradient
695 : ! (for the partial integration)
696 : ! deriv rho
697 537860 : DO idir = 1, 3
698 537860 : CALL pw_pool%create_pw(drho_r(idir, ispin))
699 : END DO
700 134465 : IF (needs_rho_g) THEN
701 132223 : IF (.NOT. ASSOCIATED(my_rho_g%pw_grid)) THEN
702 26806 : my_rho_g_local = .TRUE.
703 26806 : CALL pw_pool%create_pw(my_rho_g)
704 26806 : CALL pw_transfer(my_rho_r(ispin), my_rho_g)
705 : END IF
706 105417 : IF (.NOT. my_rho_g_local .AND. (xc_deriv_method_id == xc_deriv_spline2 .OR. &
707 : xc_deriv_method_id == xc_deriv_spline3)) THEN
708 7386 : CALL pw_pool%create_pw(my_rho_g)
709 7386 : my_rho_g_local = .TRUE.
710 7386 : CALL pw_copy(rho_g(ispin), my_rho_g)
711 : END IF
712 : END IF
713 134465 : IF (needs%laplace_rho .OR. needs%laplace_rho_spin) THEN
714 2208 : CALL pw_pool%create_pw(laplace_rho_r(ispin))
715 2208 : CALL xc_pw_laplace(my_rho_g, pw_pool, xc_deriv_method_id, laplace_rho_r(ispin), tmp_g=tmp_g)
716 : END IF
717 134465 : CALL xc_pw_gradient(my_rho_r(ispin), my_rho_g, tmp_g, drho_r(:, ispin), xc_deriv_method_id)
718 :
719 134465 : IF (needs_rho_g) THEN
720 132223 : IF (my_rho_g_local) THEN
721 34192 : my_rho_g_local = .FALSE.
722 34192 : CALL pw_pool%give_back_pw(my_rho_g)
723 : END IF
724 : END IF
725 :
726 134465 : IF (xc_deriv_method_id /= xc_deriv_pw) THEN
727 9754 : CALL pw_spline_scale_deriv(drho_r(:, ispin))
728 : END IF
729 :
730 : END IF
731 :
732 : END DO
733 :
734 191789 : IF (ASSOCIATED(tmp_g%pw_grid)) THEN
735 107963 : CALL pw_pool%give_back_pw(tmp_g)
736 : END IF
737 :
738 150106 : SELECT CASE (nspins)
739 : CASE (1)
740 150106 : IF (.NOT. do_sf) THEN
741 149998 : IF (needs%rho_1_3) THEN
742 10442 : CALL pw_pool%create_cr3d(rho_set%rho_1_3)
743 10442 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,my_rho_r)
744 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
745 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
746 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
747 : rho_set%rho_1_3(i, j, k) = MAX(my_rho_r(1)%array(i, j, k), 0.0_dp)**f13
748 : END DO
749 : END DO
750 : END DO
751 10442 : rho_set%owns%rho_1_3 = .TRUE.
752 10442 : rho_set%has%rho_1_3 = .TRUE.
753 : END IF
754 149998 : IF (needs%rho) THEN
755 149998 : rho_set%rho => my_rho_r(1)%array
756 149998 : NULLIFY (my_rho_r(1)%array)
757 149998 : rho_set%owns%rho = .TRUE.
758 149998 : rho_set%has%rho = .TRUE.
759 : END IF
760 149998 : IF (needs%norm_drho) THEN
761 84485 : CALL pw_pool%create_cr3d(rho_set%norm_drho)
762 84485 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
763 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
764 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
765 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
766 : rho_set%norm_drho(i, j, k) = SQRT( &
767 : drho_r(1, 1)%array(i, j, k)**2 + &
768 : drho_r(2, 1)%array(i, j, k)**2 + &
769 : drho_r(3, 1)%array(i, j, k)**2)
770 : END DO
771 : END DO
772 : END DO
773 84485 : rho_set%owns%norm_drho = .TRUE.
774 84485 : rho_set%has%norm_drho = .TRUE.
775 : END IF
776 149998 : IF (needs%laplace_rho) THEN
777 1104 : rho_set%laplace_rho => laplace_rho_r(1)%array
778 1104 : NULLIFY (laplace_rho_r(1)%array)
779 1104 : rho_set%owns%laplace_rho = .TRUE.
780 1104 : rho_set%has%laplace_rho = .TRUE.
781 : END IF
782 :
783 149998 : IF (needs%drho) THEN
784 324812 : DO idir = 1, 3
785 243609 : rho_set%drho(idir)%array => drho_r(idir, 1)%array
786 324812 : NULLIFY (drho_r(idir, 1)%array)
787 : END DO
788 81203 : rho_set%owns%drho = .TRUE.
789 81203 : rho_set%has%drho = .TRUE.
790 : END IF
791 : ELSE
792 108 : IF (needs%norm_drho) THEN
793 56 : CALL pw_pool%create_cr3d(rho_set%norm_drho)
794 56 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
795 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
796 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
797 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
798 : rho_set%norm_drho(i, j, k) = SQRT( &
799 : drho_r(1, 1)%array(i, j, k)**2 + &
800 : drho_r(2, 1)%array(i, j, k)**2 + &
801 : drho_r(3, 1)%array(i, j, k)**2)
802 : END DO
803 : END DO
804 : END DO
805 56 : rho_set%owns%norm_drho = .TRUE.
806 56 : rho_set%has%norm_drho = .TRUE.
807 : END IF
808 108 : IF (needs%rho_spin) THEN
809 :
810 108 : rho_set%rhoa => my_rho_r(1)%array
811 108 : NULLIFY (my_rho_r(1)%array)
812 :
813 108 : rho_set%owns%rho_spin = .TRUE.
814 108 : rho_set%has%rho_spin = .TRUE.
815 : END IF
816 108 : IF (needs%norm_drho_spin) THEN
817 56 : CALL pw_pool%create_cr3d(rho_set%norm_drhoa)
818 56 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
819 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
820 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
821 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
822 : rho_set%norm_drhoa(i, j, k) = SQRT( &
823 : drho_r(1, 1)%array(i, j, k)**2 + &
824 : drho_r(2, 1)%array(i, j, k)**2 + &
825 : drho_r(3, 1)%array(i, j, k)**2)
826 : END DO
827 : END DO
828 : END DO
829 56 : rho_set%owns%norm_drho_spin = .TRUE.
830 56 : rho_set%has%norm_drho_spin = .TRUE.
831 : END IF
832 108 : IF (needs%laplace_rho_spin) THEN
833 0 : rho_set%laplace_rhoa => laplace_rho_r(1)%array
834 0 : NULLIFY (laplace_rho_r(1)%array)
835 :
836 0 : rho_set%owns%laplace_rho_spin = .TRUE.
837 0 : rho_set%has%laplace_rho_spin = .TRUE.
838 : END IF
839 108 : IF (needs%drho_spin) THEN
840 224 : DO idir = 1, 3
841 168 : rho_set%drhoa(idir)%array => drho_r(idir, 1)%array
842 224 : NULLIFY (drho_r(idir, 1)%array)
843 : END DO
844 56 : rho_set%owns%drho_spin = .TRUE.
845 56 : rho_set%has%drho_spin = .TRUE.
846 : END IF
847 : END IF
848 : CASE (2)
849 41683 : IF (needs%rho_spin_1_3) THEN
850 1776 : CALL pw_pool%create_cr3d(rho_set%rhoa_1_3)
851 : !assume that the bounds are the same?
852 1776 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,my_rho_r)
853 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
854 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
855 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
856 : rho_set%rhoa_1_3(i, j, k) = MAX(my_rho_r(1)%array(i, j, k), 0.0_dp)**f13
857 : END DO
858 : END DO
859 : END DO
860 1776 : CALL pw_pool%create_cr3d(rho_set%rhob_1_3)
861 : !assume that the bounds are the same?
862 1776 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,my_rho_r)
863 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
864 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
865 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
866 : rho_set%rhob_1_3(i, j, k) = MAX(my_rho_r(2)%array(i, j, k), 0.0_dp)**f13
867 : END DO
868 : END DO
869 : END DO
870 1776 : rho_set%owns%rho_spin_1_3 = .TRUE.
871 1776 : rho_set%has%rho_spin_1_3 = .TRUE.
872 : END IF
873 41683 : IF (needs%norm_drho) THEN
874 :
875 23420 : CALL pw_pool%create_cr3d(rho_set%norm_drho)
876 23420 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
877 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
878 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
879 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
880 : rho_set%norm_drho(i, j, k) = SQRT( &
881 : (drho_r(1, 1)%array(i, j, k) + drho_r(1, 2)%array(i, j, k))**2 + &
882 : (drho_r(2, 1)%array(i, j, k) + drho_r(2, 2)%array(i, j, k))**2 + &
883 : (drho_r(3, 1)%array(i, j, k) + drho_r(3, 2)%array(i, j, k))**2)
884 : END DO
885 : END DO
886 : END DO
887 :
888 23420 : rho_set%owns%norm_drho = .TRUE.
889 23420 : rho_set%has%norm_drho = .TRUE.
890 : END IF
891 41683 : IF (needs%rho_spin) THEN
892 :
893 41683 : rho_set%rhoa => my_rho_r(1)%array
894 41683 : NULLIFY (my_rho_r(1)%array)
895 :
896 41683 : rho_set%rhob => my_rho_r(2)%array
897 41683 : NULLIFY (my_rho_r(2)%array)
898 :
899 41683 : rho_set%owns%rho_spin = .TRUE.
900 41683 : rho_set%has%rho_spin = .TRUE.
901 : END IF
902 41683 : IF (needs%norm_drho_spin) THEN
903 :
904 24782 : CALL pw_pool%create_cr3d(rho_set%norm_drhoa)
905 24782 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
906 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
907 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
908 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
909 : rho_set%norm_drhoa(i, j, k) = SQRT( &
910 : drho_r(1, 1)%array(i, j, k)**2 + &
911 : drho_r(2, 1)%array(i, j, k)**2 + &
912 : drho_r(3, 1)%array(i, j, k)**2)
913 : END DO
914 : END DO
915 : END DO
916 :
917 24782 : CALL pw_pool%create_cr3d(rho_set%norm_drhob)
918 24782 : rho_set%owns%norm_drho_spin = .TRUE.
919 24782 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
920 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
921 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
922 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
923 : rho_set%norm_drhob(i, j, k) = SQRT( &
924 : drho_r(1, 2)%array(i, j, k)**2 + &
925 : drho_r(2, 2)%array(i, j, k)**2 + &
926 : drho_r(3, 2)%array(i, j, k)**2)
927 : END DO
928 : END DO
929 : END DO
930 :
931 24782 : rho_set%owns%norm_drho_spin = .TRUE.
932 24782 : rho_set%has%norm_drho_spin = .TRUE.
933 : END IF
934 41683 : IF (needs%laplace_rho_spin) THEN
935 552 : rho_set%laplace_rhoa => laplace_rho_r(1)%array
936 552 : NULLIFY (laplace_rho_r(1)%array)
937 :
938 552 : rho_set%laplace_rhob => laplace_rho_r(2)%array
939 552 : NULLIFY (laplace_rho_r(2)%array)
940 :
941 552 : rho_set%owns%laplace_rho_spin = .TRUE.
942 552 : rho_set%has%laplace_rho_spin = .TRUE.
943 : END IF
944 233472 : IF (needs%drho_spin) THEN
945 86352 : DO idir = 1, 3
946 64764 : rho_set%drhoa(idir)%array => drho_r(idir, 1)%array
947 64764 : NULLIFY (drho_r(idir, 1)%array)
948 64764 : rho_set%drhob(idir)%array => drho_r(idir, 2)%array
949 86352 : NULLIFY (drho_r(idir, 2)%array)
950 : END DO
951 21588 : rho_set%owns%drho_spin = .TRUE.
952 21588 : rho_set%has%drho_spin = .TRUE.
953 : END IF
954 : END SELECT
955 : ! post cleanup
956 425261 : DO ispin = 1, nspins
957 233472 : IF (needs%laplace_rho .OR. needs%laplace_rho_spin) THEN
958 2208 : CALL pw_pool%give_back_pw(laplace_rho_r(ispin))
959 : END IF
960 1125677 : DO idir = 1, 3
961 933888 : CALL pw_pool%give_back_pw(drho_r(idir, ispin))
962 : END DO
963 : END DO
964 425261 : DO ispin = 1, nspins
965 425261 : CALL pw_pool%give_back_pw(my_rho_r(ispin))
966 : END DO
967 :
968 : ! tau part
969 191789 : IF (needs%tau .OR. needs%tau_spin) THEN
970 5796 : CPASSERT(ASSOCIATED(tau))
971 12896 : DO ispin = 1, nspins
972 198889 : CPASSERT(ASSOCIATED(tau(ispin)%array))
973 : END DO
974 : END IF
975 191789 : IF (needs%tau) THEN
976 4492 : rho_set%tau => tau(1)%array
977 4492 : rho_set%owns%tau = .FALSE.
978 4492 : rho_set%has%tau = .TRUE.
979 : END IF
980 191789 : IF (needs%tau_spin) THEN
981 1304 : rho_set%tau_a => tau(1)%array
982 1304 : rho_set%tau_b => tau(2)%array
983 1304 : rho_set%owns%tau_spin = .FALSE.
984 1304 : rho_set%has%tau_spin = .TRUE.
985 : END IF
986 :
987 191789 : CPASSERT(xc_rho_cflags_equal(rho_set%has, needs))
988 :
989 383578 : END SUBROUTINE xc_rho_set_update
990 :
991 0 : END MODULE xc_rho_set_types
|
