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 348274 : 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 348274 : IF (PRESENT(rho_cutoff)) rho_set%rho_cutoff = rho_cutoff
116 348274 : IF (PRESENT(drho_cutoff)) rho_set%drho_cutoff = drho_cutoff
117 348274 : IF (PRESENT(tau_cutoff)) rho_set%tau_cutoff = tau_cutoff
118 3482740 : rho_set%local_bounds = local_bounds
119 348274 : CALL xc_rho_cflags_setall(rho_set%owns, .TRUE.)
120 348274 : CALL xc_rho_cflags_setall(rho_set%has, .FALSE.)
121 348274 : 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 348274 : 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 348274 : IF (PRESENT(pw_pool)) THEN
135 147247 : IF (ASSOCIATED(pw_pool)) THEN
136 147247 : 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 1393096 : rho_set%local_bounds(1, :) = -HUGE(0) ! we want to crash...
143 1393096 : rho_set%local_bounds(1, :) = HUGE(0)
144 348274 : IF (rho_set%owns%rho .AND. ASSOCIATED(rho_set%rho)) THEN
145 179028 : DEALLOCATE (rho_set%rho)
146 : ELSE
147 169246 : NULLIFY (rho_set%rho)
148 : END IF
149 348274 : IF (rho_set%owns%rho_spin) THEN
150 175411 : IF (ASSOCIATED(rho_set%rhoa)) THEN
151 21871 : DEALLOCATE (rho_set%rhoa)
152 : END IF
153 175411 : IF (ASSOCIATED(rho_set%rhob)) THEN
154 21767 : DEALLOCATE (rho_set%rhob)
155 : END IF
156 : ELSE
157 172863 : NULLIFY (rho_set%rhoa, rho_set%rhob)
158 : END IF
159 348274 : IF (rho_set%owns%rho_1_3 .AND. ASSOCIATED(rho_set%rho_1_3)) THEN
160 5265 : DEALLOCATE (rho_set%rho_1_3)
161 : ELSE
162 343009 : NULLIFY (rho_set%rho_1_3)
163 : END IF
164 348274 : IF (rho_set%owns%rho_spin) THEN
165 175411 : IF (ASSOCIATED(rho_set%rhoa_1_3)) THEN
166 2452 : DEALLOCATE (rho_set%rhoa_1_3)
167 : END IF
168 175411 : IF (ASSOCIATED(rho_set%rhob_1_3)) THEN
169 2452 : DEALLOCATE (rho_set%rhob_1_3)
170 : END IF
171 : ELSE
172 172863 : NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
173 : END IF
174 348274 : IF (rho_set%owns%drho) THEN
175 744268 : DO i = 1, 3
176 744268 : IF (ASSOCIATED(rho_set%drho(i)%array)) THEN
177 315138 : DEALLOCATE (rho_set%drho(i)%array)
178 : END IF
179 : END DO
180 : ELSE
181 648828 : DO i = 1, 3
182 648828 : NULLIFY (rho_set%drho(i)%array)
183 : END DO
184 : END IF
185 348274 : IF (rho_set%owns%drho_spin) THEN
186 693772 : DO i = 1, 3
187 520329 : IF (ASSOCIATED(rho_set%drhoa(i)%array)) THEN
188 38994 : DEALLOCATE (rho_set%drhoa(i)%array)
189 : END IF
190 693772 : IF (ASSOCIATED(rho_set%drhob(i)%array)) THEN
191 38838 : DEALLOCATE (rho_set%drhob(i)%array)
192 : END IF
193 : END DO
194 : ELSE
195 699324 : DO i = 1, 3
196 699324 : NULLIFY (rho_set%drhoa(i)%array, rho_set%drhob(i)%array)
197 : END DO
198 : END IF
199 348274 : IF (rho_set%owns%laplace_rho .AND. ASSOCIATED(rho_set%laplace_rho)) THEN
200 202 : DEALLOCATE (rho_set%laplace_rho)
201 : ELSE
202 348072 : NULLIFY (rho_set%laplace_rho)
203 : END IF
204 :
205 348274 : IF (rho_set%owns%norm_drho .AND. ASSOCIATED(rho_set%norm_drho)) THEN
206 140253 : DEALLOCATE (rho_set%norm_drho)
207 : ELSE
208 208021 : NULLIFY (rho_set%norm_drho)
209 : END IF
210 348274 : IF (rho_set%owns%laplace_rho_spin) THEN
211 171463 : IF (ASSOCIATED(rho_set%laplace_rhoa)) THEN
212 50 : DEALLOCATE (rho_set%laplace_rhoa)
213 : END IF
214 171463 : IF (ASSOCIATED(rho_set%laplace_rhob)) THEN
215 50 : DEALLOCATE (rho_set%laplace_rhob)
216 : END IF
217 : ELSE
218 176811 : NULLIFY (rho_set%laplace_rhoa, rho_set%laplace_rhob)
219 : END IF
220 :
221 348274 : IF (rho_set%owns%norm_drho_spin) THEN
222 173443 : IF (ASSOCIATED(rho_set%norm_drhoa)) THEN
223 13877 : DEALLOCATE (rho_set%norm_drhoa)
224 : END IF
225 173443 : IF (ASSOCIATED(rho_set%norm_drhob)) THEN
226 13825 : DEALLOCATE (rho_set%norm_drhob)
227 : END IF
228 : ELSE
229 174831 : NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
230 : END IF
231 348274 : IF (rho_set%owns%tau .AND. ASSOCIATED(rho_set%tau)) THEN
232 2544 : DEALLOCATE (rho_set%tau)
233 : ELSE
234 345730 : NULLIFY (rho_set%tau)
235 : END IF
236 348274 : IF (rho_set%owns%tau_spin) THEN
237 171413 : IF (ASSOCIATED(rho_set%tau_a)) THEN
238 34 : DEALLOCATE (rho_set%tau_a)
239 : END IF
240 171413 : IF (ASSOCIATED(rho_set%tau_b)) THEN
241 34 : DEALLOCATE (rho_set%tau_b)
242 : END IF
243 : ELSE
244 176861 : NULLIFY (rho_set%tau_a, rho_set%tau_b)
245 : END IF
246 348274 : 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 1346730 : 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 1346730 : my_can_return_null = .FALSE.
301 1346730 : IF (PRESENT(can_return_null)) my_can_return_null = can_return_null
302 :
303 1346730 : IF (PRESENT(rho)) THEN
304 354317 : rho => rho_set%rho
305 354317 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rho))
306 : END IF
307 1346730 : IF (PRESENT(drho)) THEN
308 927396 : DO i = 1, 3
309 695547 : drho(i)%array => rho_set%drho(i)%array
310 927396 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rho_set%drho(i)%array))
311 : END DO
312 : END IF
313 1346730 : IF (PRESENT(norm_drho)) THEN
314 490323 : norm_drho => rho_set%norm_drho
315 490323 : CPASSERT(my_can_return_null .OR. ASSOCIATED(norm_drho))
316 : END IF
317 1346730 : IF (PRESENT(laplace_rho)) THEN
318 18466 : laplace_rho => rho_set%laplace_rho
319 18466 : CPASSERT(my_can_return_null .OR. ASSOCIATED(laplace_rho))
320 : END IF
321 1346730 : IF (PRESENT(rhoa)) THEN
322 175904 : rhoa => rho_set%rhoa
323 175904 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rhoa))
324 : END IF
325 1346730 : IF (PRESENT(rhob)) THEN
326 175490 : rhob => rho_set%rhob
327 175490 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rhob))
328 : END IF
329 1346730 : IF (PRESENT(drhoa)) THEN
330 746820 : DO i = 1, 3
331 560115 : drhoa(i)%array => rho_set%drhoa(i)%array
332 746820 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rho_set%drhoa(i)%array))
333 : END DO
334 : END IF
335 1346730 : IF (PRESENT(drhob)) THEN
336 745772 : DO i = 1, 3
337 559329 : drhob(i)%array => rho_set%drhob(i)%array
338 745772 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rho_set%drhob(i)%array))
339 : END DO
340 : END IF
341 1346730 : IF (PRESENT(laplace_rhoa)) THEN
342 2968 : laplace_rhoa => rho_set%laplace_rhoa
343 2968 : CPASSERT(my_can_return_null .OR. ASSOCIATED(laplace_rhoa))
344 : END IF
345 1346730 : IF (PRESENT(laplace_rhob)) THEN
346 2968 : laplace_rhob => rho_set%laplace_rhob
347 2968 : CPASSERT(my_can_return_null .OR. ASSOCIATED(laplace_rhob))
348 : END IF
349 1346730 : IF (PRESENT(norm_drhoa)) THEN
350 243139 : norm_drhoa => rho_set%norm_drhoa
351 243139 : CPASSERT(my_can_return_null .OR. ASSOCIATED(norm_drhoa))
352 : END IF
353 1346730 : IF (PRESENT(norm_drhob)) THEN
354 243139 : norm_drhob => rho_set%norm_drhob
355 243139 : CPASSERT(my_can_return_null .OR. ASSOCIATED(norm_drhob))
356 : END IF
357 1346730 : IF (PRESENT(rho_1_3)) THEN
358 23135 : rho_1_3 => rho_set%rho_1_3
359 23135 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rho_1_3))
360 : END IF
361 1346730 : IF (PRESENT(rhoa_1_3)) THEN
362 3364 : rhoa_1_3 => rho_set%rhoa_1_3
363 3364 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rhoa_1_3))
364 : END IF
365 1346730 : IF (PRESENT(rhob_1_3)) THEN
366 3364 : rhob_1_3 => rho_set%rhob_1_3
367 3364 : CPASSERT(my_can_return_null .OR. ASSOCIATED(rhob_1_3))
368 : END IF
369 1346730 : IF (PRESENT(tau)) THEN
370 21982 : tau => rho_set%tau
371 21982 : CPASSERT(my_can_return_null .OR. ASSOCIATED(tau))
372 : END IF
373 1346730 : IF (PRESENT(tau_a)) THEN
374 3200 : tau_a => rho_set%tau_a
375 3200 : CPASSERT(my_can_return_null .OR. ASSOCIATED(tau_a))
376 : END IF
377 1346730 : IF (PRESENT(tau_b)) THEN
378 3200 : tau_b => rho_set%tau_b
379 3200 : CPASSERT(my_can_return_null .OR. ASSOCIATED(tau_b))
380 : END IF
381 1346730 : IF (PRESENT(rho_cutoff)) rho_cutoff = rho_set%rho_cutoff
382 1346730 : IF (PRESENT(drho_cutoff)) drho_cutoff = rho_set%drho_cutoff
383 1346730 : IF (PRESENT(tau_cutoff)) tau_cutoff = rho_set%tau_cutoff
384 3515560 : IF (PRESENT(local_bounds)) local_bounds = rho_set%local_bounds
385 1346730 : 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 535535 : 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 535535 : $:recover(variable)
465 : #:endfor
466 :
467 107107 : END SUBROUTINE xc_rho_set_recover_pw
468 :
469 321321 : 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 321321 : IF (ASSOCIATED(rho)) THEN
477 276339 : CALL rho_pw%create(pw_grid=pw_grid, array_ptr=rho)
478 276339 : NULLIFY (rho)
479 44982 : ELSE IF (PRESENT(rhoa) .AND. PRESENT(rhob)) THEN
480 44982 : CALL pw_pool%create_pw(rho_pw)
481 44982 : !$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 321321 : 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 332752 : 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 332752 : IF (rho_set%owns%rho) THEN
505 300667 : CALL pw_pool%give_back_cr3d(rho_set%rho)
506 : ELSE
507 32085 : NULLIFY (rho_set%rho)
508 : END IF
509 332752 : IF (rho_set%owns%rho_1_3) THEN
510 186577 : CALL pw_pool%give_back_cr3d(rho_set%rho_1_3)
511 : ELSE
512 146175 : NULLIFY (rho_set%rho_1_3)
513 : END IF
514 332752 : IF (rho_set%owns%drho) THEN
515 973936 : DO idir = 1, 3
516 973936 : CALL pw_pool%give_back_cr3d(rho_set%drho(idir)%array)
517 : END DO
518 : ELSE
519 357072 : DO idir = 1, 3
520 357072 : NULLIFY (rho_set%drho(idir)%array)
521 : END DO
522 : END IF
523 332752 : IF (rho_set%owns%norm_drho) THEN
524 265772 : CALL pw_pool%give_back_cr3d(rho_set%norm_drho)
525 : ELSE
526 66980 : NULLIFY (rho_set%norm_drho)
527 : END IF
528 332752 : IF (rho_set%owns%laplace_rho) THEN
529 177637 : CALL pw_pool%give_back_cr3d(rho_set%laplace_rho)
530 : ELSE
531 155115 : NULLIFY (rho_set%laplace_rho)
532 : END IF
533 332752 : IF (rho_set%owns%tau) THEN
534 176861 : CALL pw_pool%give_back_cr3d(rho_set%tau)
535 : ELSE
536 155891 : NULLIFY (rho_set%tau)
537 : END IF
538 332752 : IF (rho_set%owns%rho_spin) THEN
539 208946 : CALL pw_pool%give_back_cr3d(rho_set%rhoa)
540 208946 : CALL pw_pool%give_back_cr3d(rho_set%rhob)
541 : ELSE
542 123806 : NULLIFY (rho_set%rhoa, rho_set%rhob)
543 : END IF
544 332752 : IF (rho_set%owns%rho_spin_1_3) THEN
545 178621 : CALL pw_pool%give_back_cr3d(rho_set%rhoa_1_3)
546 178621 : CALL pw_pool%give_back_cr3d(rho_set%rhob_1_3)
547 : ELSE
548 154131 : NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
549 : END IF
550 332752 : IF (rho_set%owns%drho_spin) THEN
551 776956 : DO idir = 1, 3
552 582717 : CALL pw_pool%give_back_cr3d(rho_set%drhoa(idir)%array)
553 776956 : CALL pw_pool%give_back_cr3d(rho_set%drhob(idir)%array)
554 : END DO
555 : ELSE
556 554052 : DO idir = 1, 3
557 554052 : NULLIFY (rho_set%drhoa(idir)%array, rho_set%drhob(idir)%array)
558 : END DO
559 : END IF
560 332752 : IF (rho_set%owns%laplace_rho_spin) THEN
561 177161 : CALL pw_pool%give_back_cr3d(rho_set%laplace_rhoa)
562 177161 : CALL pw_pool%give_back_cr3d(rho_set%laplace_rhob)
563 : ELSE
564 155591 : NULLIFY (rho_set%laplace_rhoa, rho_set%laplace_rhob)
565 : END IF
566 332752 : IF (rho_set%owns%norm_drho_spin) THEN
567 197267 : CALL pw_pool%give_back_cr3d(rho_set%norm_drhoa)
568 197267 : CALL pw_pool%give_back_cr3d(rho_set%norm_drhob)
569 : ELSE
570 135485 : NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
571 : END IF
572 332752 : IF (rho_set%owns%tau_spin) THEN
573 176861 : CALL pw_pool%give_back_cr3d(rho_set%tau_a)
574 176861 : CALL pw_pool%give_back_cr3d(rho_set%tau_b)
575 : ELSE
576 155891 : NULLIFY (rho_set%tau_a, rho_set%tau_b)
577 : END IF
578 :
579 332752 : CALL xc_rho_cflags_setall(rho_set%has, .FALSE.)
580 332752 : CALL xc_rho_cflags_setall(rho_set%owns, .FALSE.)
581 :
582 332752 : 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 185505 : 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 556515 : TYPE(pw_r3d_rs_type), DIMENSION(2) :: laplace_rho_r
616 1669545 : TYPE(pw_r3d_rs_type), DIMENSION(3, 2) :: drho_r
617 : TYPE(pw_c1d_gs_type) :: my_rho_g, tmp_g
618 556515 : TYPE(pw_r3d_rs_type), DIMENSION(2) :: my_rho_r
619 :
620 185505 : do_sf = .FALSE.
621 15712 : IF (PRESENT(spinflip)) do_sf = spinflip
622 :
623 1855050 : 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 185505 : nspins = SIZE(rho_r)
627 1855050 : rho_set%local_bounds = rho_r(1)%pw_grid%bounds_local
628 185505 : rho_cutoff = 0.5*rho_set%rho_cutoff
629 :
630 185505 : my_rho_g_local = .FALSE.
631 : ! some checks
632 149526 : SELECT CASE (nspins)
633 : CASE (1)
634 149526 : IF (.NOT. do_sf) THEN
635 149422 : CPASSERT(.NOT. needs%rho_spin)
636 149422 : CPASSERT(.NOT. needs%drho_spin)
637 149422 : CPASSERT(.NOT. needs%norm_drho_spin)
638 149422 : CPASSERT(.NOT. needs%rho_spin_1_3)
639 149422 : CPASSERT(.NOT. needs%tau_spin)
640 149422 : CPASSERT(.NOT. needs%laplace_rho_spin)
641 : ELSE
642 104 : CPASSERT(.NOT. needs%rho)
643 104 : CPASSERT(.NOT. needs%drho)
644 104 : CPASSERT(.NOT. needs%rho_1_3)
645 104 : CPASSERT(.NOT. needs%tau)
646 104 : CPASSERT(.NOT. needs%laplace_rho)
647 : END IF
648 : CASE (2)
649 35979 : CPASSERT(.NOT. needs%rho)
650 35979 : CPASSERT(.NOT. needs%drho)
651 35979 : CPASSERT(.NOT. needs%rho_1_3)
652 35979 : CPASSERT(.NOT. needs%tau)
653 35979 : CPASSERT(.NOT. needs%laplace_rho)
654 : CASE default
655 185505 : CPABORT("Unknown number of spin states")
656 : END SELECT
657 :
658 185505 : CALL xc_rho_set_clean(rho_set, pw_pool=pw_pool)
659 :
660 185505 : 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 185505 : needs_laplace)
664 : needs_rho_g = ((xc_deriv_method_id == xc_deriv_spline3 .OR. &
665 : xc_deriv_method_id == xc_deriv_spline2 .OR. &
666 185505 : xc_deriv_method_id == xc_deriv_pw)) .AND. (gradient_f .OR. needs_laplace)
667 107479 : 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 185505 : IF (needs_rho_g) THEN
675 105363 : CALL pw_pool%create_pw(tmp_g)
676 : END IF
677 406989 : DO ispin = 1, nspins
678 221484 : CALL pw_pool%create_pw(my_rho_r(ispin))
679 : ! introduce a smoothing kernel on the density
680 221484 : IF (xc_rho_smooth_id == xc_rho_no_smooth) THEN
681 220908 : IF (needs_rho_g) THEN
682 127147 : IF (ASSOCIATED(rho_g)) THEN
683 109163 : my_rho_g_local = .FALSE.
684 109163 : my_rho_g = rho_g(ispin)
685 : END IF
686 : END IF
687 :
688 220908 : 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 406989 : 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 518420 : DO idir = 1, 3
698 518420 : CALL pw_pool%create_pw(drho_r(idir, ispin))
699 : END DO
700 129605 : IF (needs_rho_g) THEN
701 127225 : IF (.NOT. ASSOCIATED(my_rho_g%pw_grid)) THEN
702 18062 : my_rho_g_local = .TRUE.
703 18062 : CALL pw_pool%create_pw(my_rho_g)
704 18062 : CALL pw_transfer(my_rho_r(ispin), my_rho_g)
705 : END IF
706 109163 : 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 7514 : CALL pw_pool%create_pw(my_rho_g)
709 7514 : my_rho_g_local = .TRUE.
710 7514 : CALL pw_copy(rho_g(ispin), my_rho_g)
711 : END IF
712 : END IF
713 129605 : IF (needs%laplace_rho .OR. needs%laplace_rho_spin) THEN
714 1678 : CALL pw_pool%create_pw(laplace_rho_r(ispin))
715 1678 : 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 129605 : CALL xc_pw_gradient(my_rho_r(ispin), my_rho_g, tmp_g, drho_r(:, ispin), xc_deriv_method_id)
718 :
719 129605 : IF (needs_rho_g) THEN
720 127225 : IF (my_rho_g_local) THEN
721 25576 : my_rho_g_local = .FALSE.
722 25576 : CALL pw_pool%give_back_pw(my_rho_g)
723 : END IF
724 : END IF
725 :
726 129605 : IF (xc_deriv_method_id /= xc_deriv_pw) THEN
727 9960 : CALL pw_spline_scale_deriv(drho_r(:, ispin))
728 : END IF
729 :
730 : END IF
731 :
732 : END DO
733 :
734 185505 : IF (ASSOCIATED(tmp_g%pw_grid)) THEN
735 105363 : CALL pw_pool%give_back_pw(tmp_g)
736 : END IF
737 :
738 149526 : SELECT CASE (nspins)
739 : CASE (1)
740 149526 : IF (.NOT. do_sf) THEN
741 149422 : IF (needs%rho_1_3) THEN
742 10950 : CALL pw_pool%create_cr3d(rho_set%rho_1_3)
743 10950 : !$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 10950 : rho_set%owns%rho_1_3 = .TRUE.
752 10950 : rho_set%has%rho_1_3 = .TRUE.
753 : END IF
754 149422 : IF (needs%rho) THEN
755 149422 : rho_set%rho => my_rho_r(1)%array
756 149422 : NULLIFY (my_rho_r(1)%array)
757 149422 : rho_set%owns%rho = .TRUE.
758 149422 : rho_set%has%rho = .TRUE.
759 : END IF
760 149422 : IF (needs%norm_drho) THEN
761 84447 : CALL pw_pool%create_cr3d(rho_set%norm_drho)
762 84447 : !$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 84447 : rho_set%owns%norm_drho = .TRUE.
774 84447 : rho_set%has%norm_drho = .TRUE.
775 : END IF
776 149422 : IF (needs%laplace_rho) THEN
777 978 : rho_set%laplace_rho => laplace_rho_r(1)%array
778 978 : NULLIFY (laplace_rho_r(1)%array)
779 978 : rho_set%owns%laplace_rho = .TRUE.
780 978 : rho_set%has%laplace_rho = .TRUE.
781 : END IF
782 :
783 149422 : IF (needs%drho) THEN
784 325108 : DO idir = 1, 3
785 243831 : rho_set%drho(idir)%array => drho_r(idir, 1)%array
786 325108 : NULLIFY (drho_r(idir, 1)%array)
787 : END DO
788 81277 : rho_set%owns%drho = .TRUE.
789 81277 : rho_set%has%drho = .TRUE.
790 : END IF
791 : ELSE
792 104 : IF (needs%norm_drho) THEN
793 52 : CALL pw_pool%create_cr3d(rho_set%norm_drho)
794 52 : !$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 52 : rho_set%owns%norm_drho = .TRUE.
806 52 : rho_set%has%norm_drho = .TRUE.
807 : END IF
808 104 : IF (needs%rho_spin) THEN
809 :
810 104 : rho_set%rhoa => my_rho_r(1)%array
811 104 : NULLIFY (my_rho_r(1)%array)
812 :
813 104 : rho_set%owns%rho_spin = .TRUE.
814 104 : rho_set%has%rho_spin = .TRUE.
815 : END IF
816 104 : IF (needs%norm_drho_spin) THEN
817 52 : CALL pw_pool%create_cr3d(rho_set%norm_drhoa)
818 52 : !$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 52 : rho_set%owns%norm_drho_spin = .TRUE.
830 52 : rho_set%has%norm_drho_spin = .TRUE.
831 : END IF
832 104 : 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 104 : IF (needs%drho_spin) THEN
840 208 : DO idir = 1, 3
841 156 : rho_set%drhoa(idir)%array => drho_r(idir, 1)%array
842 208 : NULLIFY (drho_r(idir, 1)%array)
843 : END DO
844 52 : rho_set%owns%drho_spin = .TRUE.
845 52 : rho_set%has%drho_spin = .TRUE.
846 : END IF
847 : END IF
848 : CASE (2)
849 35979 : IF (needs%rho_spin_1_3) THEN
850 1772 : CALL pw_pool%create_cr3d(rho_set%rhoa_1_3)
851 : !assume that the bounds are the same?
852 1772 : !$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 1772 : CALL pw_pool%create_cr3d(rho_set%rhob_1_3)
861 : !assume that the bounds are the same?
862 1772 : !$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 1772 : rho_set%owns%rho_spin_1_3 = .TRUE.
871 1772 : rho_set%has%rho_spin_1_3 = .TRUE.
872 : END IF
873 35979 : IF (needs%norm_drho) THEN
874 :
875 21076 : CALL pw_pool%create_cr3d(rho_set%norm_drho)
876 21076 : !$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 21076 : rho_set%owns%norm_drho = .TRUE.
889 21076 : rho_set%has%norm_drho = .TRUE.
890 : END IF
891 35979 : IF (needs%rho_spin) THEN
892 :
893 35979 : rho_set%rhoa => my_rho_r(1)%array
894 35979 : NULLIFY (my_rho_r(1)%array)
895 :
896 35979 : rho_set%rhob => my_rho_r(2)%array
897 35979 : NULLIFY (my_rho_r(2)%array)
898 :
899 35979 : rho_set%owns%rho_spin = .TRUE.
900 35979 : rho_set%has%rho_spin = .TRUE.
901 : END IF
902 35979 : IF (needs%norm_drho_spin) THEN
903 :
904 22384 : CALL pw_pool%create_cr3d(rho_set%norm_drhoa)
905 22384 : !$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 22384 : CALL pw_pool%create_cr3d(rho_set%norm_drhob)
918 22384 : rho_set%owns%norm_drho_spin = .TRUE.
919 22384 : !$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 22384 : rho_set%owns%norm_drho_spin = .TRUE.
932 22384 : rho_set%has%norm_drho_spin = .TRUE.
933 : END IF
934 35979 : IF (needs%laplace_rho_spin) THEN
935 350 : rho_set%laplace_rhoa => laplace_rho_r(1)%array
936 350 : NULLIFY (laplace_rho_r(1)%array)
937 :
938 350 : rho_set%laplace_rhob => laplace_rho_r(2)%array
939 350 : NULLIFY (laplace_rho_r(2)%array)
940 :
941 350 : rho_set%owns%laplace_rho_spin = .TRUE.
942 350 : rho_set%has%laplace_rho_spin = .TRUE.
943 : END IF
944 221484 : IF (needs%drho_spin) THEN
945 77424 : DO idir = 1, 3
946 58068 : rho_set%drhoa(idir)%array => drho_r(idir, 1)%array
947 58068 : NULLIFY (drho_r(idir, 1)%array)
948 58068 : rho_set%drhob(idir)%array => drho_r(idir, 2)%array
949 77424 : NULLIFY (drho_r(idir, 2)%array)
950 : END DO
951 19356 : rho_set%owns%drho_spin = .TRUE.
952 19356 : rho_set%has%drho_spin = .TRUE.
953 : END IF
954 : END SELECT
955 : ! post cleanup
956 406989 : DO ispin = 1, nspins
957 221484 : IF (needs%laplace_rho .OR. needs%laplace_rho_spin) THEN
958 1678 : CALL pw_pool%give_back_pw(laplace_rho_r(ispin))
959 : END IF
960 1071441 : DO idir = 1, 3
961 885936 : CALL pw_pool%give_back_pw(drho_r(idir, ispin))
962 : END DO
963 : END DO
964 406989 : DO ispin = 1, nspins
965 406989 : CALL pw_pool%give_back_pw(my_rho_r(ispin))
966 : END DO
967 :
968 : ! tau part
969 185505 : IF (needs%tau .OR. needs%tau_spin) THEN
970 4848 : CPASSERT(ASSOCIATED(tau))
971 10566 : DO ispin = 1, nspins
972 191223 : CPASSERT(ASSOCIATED(tau(ispin)%array))
973 : END DO
974 : END IF
975 185505 : IF (needs%tau) THEN
976 3978 : rho_set%tau => tau(1)%array
977 3978 : rho_set%owns%tau = .FALSE.
978 3978 : rho_set%has%tau = .TRUE.
979 : END IF
980 185505 : IF (needs%tau_spin) THEN
981 870 : rho_set%tau_a => tau(1)%array
982 870 : rho_set%tau_b => tau(2)%array
983 870 : rho_set%owns%tau_spin = .FALSE.
984 870 : rho_set%has%tau_spin = .TRUE.
985 : END IF
986 :
987 185505 : CPASSERT(xc_rho_cflags_equal(rho_set%has, needs))
988 :
989 371010 : END SUBROUTINE xc_rho_set_update
990 :
991 0 : END MODULE xc_rho_set_types
|
