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 This module defines the grid data type and some basic operations on it
10 : !> \note
11 : !> pw_grid_create : set the defaults
12 : !> pw_grid_release : release all memory connected to type
13 : !> pw_grid_setup : main routine to set up a grid
14 : !> input: cell (the box for the grid)
15 : !> pw_grid (the grid; pw_grid%grid_span has to be set)
16 : !> cutoff (optional, if not given pw_grid%bounds has to be set)
17 : !> pe_group (optional, if not given we have a local grid)
18 : !>
19 : !> if no cutoff or a negative cutoff is given, all g-vectors
20 : !> in the box are included (no spherical cutoff)
21 : !>
22 : !> for a distributed setup the array in para rs_dims has to
23 : !> be initialized
24 : !> output: pw_grid
25 : !>
26 : !> pw_grid_change : updates g-vectors after a change of the box
27 : !> \par History
28 : !> JGH (20-12-2000) : Adapted for parallel use
29 : !> JGH (07-02-2001) : Added constructor and destructor routines
30 : !> JGH (21-02-2003) : Generalized reference grid concept
31 : !> JGH (19-11-2007) : Refactoring and modularization
32 : !> JGH (21-12-2007) : pw_grid_setup refactoring
33 : !> \author apsi
34 : !> CJM
35 : ! **************************************************************************************************
36 : MODULE pw_grids
37 : USE ISO_C_BINDING, ONLY: C_F_POINTER,&
38 : C_LOC,&
39 : C_PTR,&
40 : C_SIZE_T
41 : USE kinds, ONLY: dp,&
42 : int_8,&
43 : int_size
44 : USE mathconstants, ONLY: twopi
45 : USE mathlib, ONLY: det_3x3,&
46 : inv_3x3
47 : USE message_passing, ONLY: mp_comm_self,&
48 : mp_comm_type,&
49 : mp_dims_create
50 : USE offload_api, ONLY: offload_activate_chosen_device,&
51 : offload_free_pinned_mem,&
52 : offload_malloc_pinned_mem
53 : USE pw_grid_info, ONLY: pw_find_cutoff,&
54 : pw_grid_bounds_from_n,&
55 : pw_grid_init_setup
56 : USE pw_grid_types, ONLY: FULLSPACE,&
57 : HALFSPACE,&
58 : PW_MODE_DISTRIBUTED,&
59 : PW_MODE_LOCAL,&
60 : map_pn,&
61 : pw_grid_type
62 : USE util, ONLY: get_limit,&
63 : sort
64 : #include "../base/base_uses.f90"
65 :
66 : IMPLICIT NONE
67 :
68 : PRIVATE
69 : PUBLIC :: pw_grid_create, pw_grid_retain, pw_grid_release
70 : PUBLIC :: get_pw_grid_info, pw_grid_compare
71 : PUBLIC :: pw_grid_change
72 :
73 : INTEGER :: grid_tag = 0
74 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_grids'
75 :
76 : ! Distribution in g-space can be
77 : INTEGER, PARAMETER, PUBLIC :: do_pw_grid_blocked_false = 0, &
78 : do_pw_grid_blocked_true = 1, &
79 : do_pw_grid_blocked_free = 2
80 :
81 : INTERFACE pw_grid_create
82 : MODULE PROCEDURE pw_grid_create_local
83 : MODULE PROCEDURE pw_grid_create_extended
84 : END INTERFACE
85 :
86 : CONTAINS
87 :
88 : ! **************************************************************************************************
89 : !> \brief Initialize a PW grid with bounds only (used by some routines)
90 : !> \param pw_grid ...
91 : !> \param bounds ...
92 : !> \par History
93 : !> JGH (21-Feb-2003) : initialize pw_grid%reference
94 : !> \author JGH (7-Feb-2001) & fawzi
95 : ! **************************************************************************************************
96 92762 : SUBROUTINE pw_grid_create_local(pw_grid, bounds)
97 :
98 : TYPE(pw_grid_type), POINTER :: pw_grid
99 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bounds
100 :
101 : INTEGER, DIMENSION(2) :: rs_dims
102 :
103 92762 : CPASSERT(.NOT. ASSOCIATED(pw_grid))
104 5472958 : ALLOCATE (pw_grid)
105 927620 : pw_grid%bounds = bounds
106 927620 : pw_grid%bounds_local = bounds
107 371048 : pw_grid%npts = bounds(2, :) - bounds(1, :) + 1
108 371048 : pw_grid%npts_local = pw_grid%npts
109 371048 : pw_grid%ngpts = PRODUCT(INT(pw_grid%npts, KIND=int_8))
110 92762 : pw_grid%ngpts_cut = pw_grid%ngpts
111 371048 : pw_grid%ngpts_local = PRODUCT(pw_grid%npts)
112 92762 : pw_grid%ngpts_cut_local = pw_grid%ngpts_local
113 92762 : pw_grid%grid_span = FULLSPACE
114 92762 : pw_grid%para%mode = PW_MODE_LOCAL
115 92762 : pw_grid%reference = 0
116 92762 : pw_grid%ref_count = 1
117 92762 : NULLIFY (pw_grid%g)
118 92762 : NULLIFY (pw_grid%gsq)
119 92762 : NULLIFY (pw_grid%g_hatmap)
120 92762 : NULLIFY (pw_grid%gidx)
121 92762 : NULLIFY (pw_grid%grays)
122 :
123 : ! assign a unique tag to this grid
124 92762 : grid_tag = grid_tag + 1
125 92762 : pw_grid%id_nr = grid_tag
126 :
127 : ! parallel info
128 278286 : rs_dims = 1
129 92762 : CALL pw_grid%para%group%create(mp_comm_self, 2, rs_dims)
130 92762 : IF (pw_grid%para%group%num_pe > 1) THEN
131 0 : pw_grid%para%mode = PW_MODE_DISTRIBUTED
132 : ELSE
133 92762 : pw_grid%para%mode = PW_MODE_LOCAL
134 : END IF
135 :
136 92762 : END SUBROUTINE pw_grid_create_local
137 :
138 : ! **************************************************************************************************
139 : !> \brief Check if two pw_grids are equal
140 : !> \param grida ...
141 : !> \param gridb ...
142 : !> \return ...
143 : !> \par History
144 : !> none
145 : !> \author JGH (14-Feb-2001)
146 : ! **************************************************************************************************
147 9388880 : FUNCTION pw_grid_compare(grida, gridb) RESULT(equal)
148 :
149 : TYPE(pw_grid_type), INTENT(IN) :: grida, gridb
150 : LOGICAL :: equal
151 :
152 : !------------------------------------------------------------------------------
153 :
154 9388880 : IF (grida%id_nr == gridb%id_nr) THEN
155 : equal = .TRUE.
156 : ELSE
157 : ! for the moment all grids with different identifiers are considered as not equal
158 : ! later we can get this more relaxed
159 1368 : equal = .FALSE.
160 : END IF
161 :
162 9388880 : END FUNCTION pw_grid_compare
163 :
164 : ! **************************************************************************************************
165 : !> \brief Access to information stored in the pw_grid_type
166 : !> \param pw_grid ...
167 : !> \param id_nr ...
168 : !> \param mode ...
169 : !> \param vol ...
170 : !> \param dvol ...
171 : !> \param npts ...
172 : !> \param ngpts ...
173 : !> \param ngpts_cut ...
174 : !> \param dr ...
175 : !> \param cutoff ...
176 : !> \param orthorhombic ...
177 : !> \param gvectors ...
178 : !> \param gsquare ...
179 : !> \par History
180 : !> none
181 : !> \author JGH (17-Nov-2007)
182 : ! **************************************************************************************************
183 78730 : SUBROUTINE get_pw_grid_info(pw_grid, id_nr, mode, vol, dvol, npts, ngpts, &
184 : ngpts_cut, dr, cutoff, orthorhombic, gvectors, gsquare)
185 :
186 : TYPE(pw_grid_type), INTENT(IN) :: pw_grid
187 : INTEGER, INTENT(OUT), OPTIONAL :: id_nr, mode
188 : REAL(dp), INTENT(OUT), OPTIONAL :: vol, dvol
189 : INTEGER, DIMENSION(3), INTENT(OUT), OPTIONAL :: npts
190 : INTEGER(int_8), INTENT(OUT), OPTIONAL :: ngpts, ngpts_cut
191 : REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: dr
192 : REAL(dp), INTENT(OUT), OPTIONAL :: cutoff
193 : LOGICAL, INTENT(OUT), OPTIONAL :: orthorhombic
194 : REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: gvectors
195 : REAL(dp), DIMENSION(:), OPTIONAL, POINTER :: gsquare
196 :
197 78730 : CPASSERT(pw_grid%ref_count > 0)
198 :
199 78730 : IF (PRESENT(id_nr)) id_nr = pw_grid%id_nr
200 78730 : IF (PRESENT(mode)) mode = pw_grid%para%mode
201 78730 : IF (PRESENT(vol)) vol = pw_grid%vol
202 78730 : IF (PRESENT(dvol)) dvol = pw_grid%dvol
203 393362 : IF (PRESENT(npts)) npts(1:3) = pw_grid%npts(1:3)
204 78730 : IF (PRESENT(ngpts)) ngpts = pw_grid%ngpts
205 78730 : IF (PRESENT(ngpts_cut)) ngpts_cut = pw_grid%ngpts_cut
206 78882 : IF (PRESENT(dr)) dr = pw_grid%dr
207 78730 : IF (PRESENT(cutoff)) cutoff = pw_grid%cutoff
208 78730 : IF (PRESENT(orthorhombic)) orthorhombic = pw_grid%orthorhombic
209 78730 : IF (PRESENT(gvectors)) gvectors => pw_grid%g
210 78730 : IF (PRESENT(gsquare)) gsquare => pw_grid%gsq
211 :
212 78730 : END SUBROUTINE get_pw_grid_info
213 :
214 : ! **************************************************************************************************
215 : !> \brief Set some information stored in the pw_grid_type
216 : !> \param pw_grid ...
217 : !> \param grid_span ...
218 : !> \param npts ...
219 : !> \param bounds ...
220 : !> \param cutoff ...
221 : !> \param spherical ...
222 : !> \par History
223 : !> none
224 : !> \author JGH (19-Nov-2007)
225 : ! **************************************************************************************************
226 87540 : SUBROUTINE set_pw_grid_info(pw_grid, grid_span, npts, bounds, cutoff, spherical)
227 :
228 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
229 : INTEGER, INTENT(in), OPTIONAL :: grid_span
230 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: npts
231 : INTEGER, DIMENSION(2, 3), INTENT(IN), OPTIONAL :: bounds
232 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: cutoff
233 : LOGICAL, INTENT(IN), OPTIONAL :: spherical
234 :
235 87540 : CPASSERT(pw_grid%ref_count > 0)
236 :
237 87540 : IF (PRESENT(grid_span)) THEN
238 42946 : pw_grid%grid_span = grid_span
239 : END IF
240 87540 : IF (PRESENT(bounds) .AND. PRESENT(npts)) THEN
241 0 : pw_grid%bounds = bounds
242 0 : pw_grid%npts = npts
243 0 : CPASSERT(ALL(npts == bounds(2, :) - bounds(1, :) + 1))
244 87540 : ELSE IF (PRESENT(bounds)) THEN
245 11750 : pw_grid%bounds = bounds
246 4700 : pw_grid%npts = bounds(2, :) - bounds(1, :) + 1
247 86365 : ELSE IF (PRESENT(npts)) THEN
248 173676 : pw_grid%npts = npts
249 434190 : pw_grid%bounds = pw_grid_bounds_from_n(npts)
250 : END IF
251 87540 : IF (PRESENT(cutoff)) THEN
252 44594 : pw_grid%cutoff = cutoff
253 44594 : IF (PRESENT(spherical)) THEN
254 44594 : pw_grid%spherical = spherical
255 : ELSE
256 0 : pw_grid%spherical = .FALSE.
257 : END IF
258 : END IF
259 :
260 87540 : END SUBROUTINE set_pw_grid_info
261 :
262 : ! **************************************************************************************************
263 : !> \brief sets up a pw_grid
264 : !> \param pw_grid ...
265 : !> \param mp_comm ...
266 : !> \param cell_hmat ...
267 : !> \param grid_span ...
268 : !> \param cutoff ...
269 : !> \param bounds ...
270 : !> \param bounds_local ...
271 : !> \param npts ...
272 : !> \param spherical ...
273 : !> \param odd ...
274 : !> \param fft_usage ...
275 : !> \param ncommensurate ...
276 : !> \param icommensurate ...
277 : !> \param blocked ...
278 : !> \param ref_grid ...
279 : !> \param rs_dims ...
280 : !> \param iounit ...
281 : !> \author JGH (21-Dec-2007)
282 : !> \note
283 : !> this is the function that should be used in the future
284 : ! **************************************************************************************************
285 44594 : SUBROUTINE pw_grid_create_extended(pw_grid, mp_comm, cell_hmat, grid_span, cutoff, bounds, bounds_local, npts, &
286 : spherical, odd, fft_usage, ncommensurate, icommensurate, blocked, ref_grid, &
287 : rs_dims, iounit)
288 :
289 : TYPE(pw_grid_type), POINTER :: pw_grid
290 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: cell_hmat
291 :
292 : CLASS(mp_comm_type), INTENT(IN) :: mp_comm
293 : INTEGER, INTENT(in), OPTIONAL :: grid_span
294 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: cutoff
295 : INTEGER, DIMENSION(2, 3), INTENT(IN), OPTIONAL :: bounds, bounds_local
296 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: npts
297 : LOGICAL, INTENT(in), OPTIONAL :: spherical, odd, fft_usage
298 : INTEGER, INTENT(in), OPTIONAL :: ncommensurate, icommensurate, blocked
299 : TYPE(pw_grid_type), INTENT(IN), OPTIONAL :: ref_grid
300 : INTEGER, DIMENSION(2), INTENT(in), OPTIONAL :: rs_dims
301 : INTEGER, INTENT(in), OPTIONAL :: iounit
302 :
303 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_setup'
304 :
305 : INTEGER :: handle, my_icommensurate, &
306 : my_ncommensurate
307 : INTEGER, DIMENSION(3) :: n
308 : LOGICAL :: my_fft_usage, my_odd, my_spherical
309 : REAL(KIND=dp) :: cell_deth, my_cutoff
310 : REAL(KIND=dp), DIMENSION(3, 3) :: cell_h_inv
311 :
312 44594 : CALL timeset(routineN, handle)
313 :
314 44594 : CPASSERT(.NOT. ASSOCIATED(pw_grid))
315 2318888 : ALLOCATE (pw_grid)
316 445940 : pw_grid%bounds = 0
317 44594 : pw_grid%cutoff = 0.0_dp
318 44594 : pw_grid%grid_span = FULLSPACE
319 : pw_grid%para%mode = PW_MODE_LOCAL
320 : pw_grid%reference = 0
321 44594 : pw_grid%ref_count = 1
322 44594 : NULLIFY (pw_grid%g)
323 44594 : NULLIFY (pw_grid%gsq)
324 44594 : NULLIFY (pw_grid%g_hatmap)
325 44594 : NULLIFY (pw_grid%gidx)
326 44594 : NULLIFY (pw_grid%grays)
327 :
328 : ! assign a unique tag to this grid
329 44594 : grid_tag = grid_tag + 1
330 44594 : pw_grid%id_nr = grid_tag
331 :
332 : ! parallel info
333 44594 : IF (mp_comm%num_pe > 1) THEN
334 39291 : pw_grid%para%mode = PW_MODE_DISTRIBUTED
335 : ELSE
336 : pw_grid%para%mode = PW_MODE_LOCAL
337 : END IF
338 :
339 44594 : cell_deth = ABS(det_3x3(cell_hmat))
340 44594 : IF (cell_deth < 1.0E-10_dp) THEN
341 : CALL cp_abort(__LOCATION__, &
342 : "An invalid set of cell vectors was specified. "// &
343 0 : "The determinant det(h) is too small")
344 : END IF
345 44594 : cell_h_inv = inv_3x3(cell_hmat)
346 :
347 44594 : IF (PRESENT(grid_span)) THEN
348 42946 : CALL set_pw_grid_info(pw_grid, grid_span=grid_span)
349 : END IF
350 :
351 44594 : IF (PRESENT(spherical)) THEN
352 44332 : my_spherical = spherical
353 : ELSE
354 262 : my_spherical = .FALSE.
355 : END IF
356 :
357 44594 : IF (PRESENT(odd)) THEN
358 39235 : my_odd = odd
359 : ELSE
360 5359 : my_odd = .FALSE.
361 : END IF
362 :
363 44594 : IF (PRESENT(fft_usage)) THEN
364 44428 : my_fft_usage = fft_usage
365 : ELSE
366 166 : my_fft_usage = .FALSE.
367 : END IF
368 :
369 44594 : IF (PRESENT(ncommensurate)) THEN
370 39171 : my_ncommensurate = ncommensurate
371 39171 : IF (PRESENT(icommensurate)) THEN
372 39171 : my_icommensurate = icommensurate
373 : ELSE
374 0 : my_icommensurate = MIN(1, ncommensurate)
375 : END IF
376 : ELSE
377 5423 : my_ncommensurate = 0
378 5423 : my_icommensurate = 1
379 : END IF
380 :
381 44594 : IF (PRESENT(bounds)) THEN
382 1175 : IF (PRESENT(cutoff)) THEN
383 : CALL set_pw_grid_info(pw_grid, bounds=bounds, cutoff=cutoff, &
384 180 : spherical=my_spherical)
385 : ELSE
386 3980 : n = bounds(2, :) - bounds(1, :) + 1
387 995 : my_cutoff = pw_find_cutoff(n, cell_h_inv)
388 995 : my_cutoff = 0.5_dp*my_cutoff*my_cutoff
389 : CALL set_pw_grid_info(pw_grid, bounds=bounds, cutoff=my_cutoff, &
390 995 : spherical=my_spherical)
391 : END IF
392 43419 : ELSE IF (PRESENT(npts)) THEN
393 2646 : n = npts
394 2646 : IF (PRESENT(cutoff)) THEN
395 24 : my_cutoff = cutoff
396 : ELSE
397 2622 : my_cutoff = pw_find_cutoff(npts, cell_h_inv)
398 2622 : my_cutoff = 0.5_dp*my_cutoff*my_cutoff
399 : END IF
400 2646 : IF (my_fft_usage) THEN
401 : n = pw_grid_init_setup(cell_hmat, cutoff=my_cutoff, &
402 : spherical=my_spherical, odd=my_odd, fft_usage=my_fft_usage, &
403 : ncommensurate=my_ncommensurate, icommensurate=my_icommensurate, &
404 10568 : ref_grid=ref_grid, n_orig=n)
405 : END IF
406 : CALL set_pw_grid_info(pw_grid, npts=n, cutoff=my_cutoff, &
407 2646 : spherical=my_spherical)
408 40773 : ELSE IF (PRESENT(cutoff)) THEN
409 : n = pw_grid_init_setup(cell_hmat, cutoff=cutoff, &
410 : spherical=my_spherical, odd=my_odd, fft_usage=my_fft_usage, &
411 : ncommensurate=my_ncommensurate, icommensurate=my_icommensurate, &
412 40773 : ref_grid=ref_grid)
413 : CALL set_pw_grid_info(pw_grid, npts=n, cutoff=cutoff, &
414 40773 : spherical=my_spherical)
415 : ELSE
416 0 : CPABORT("BOUNDS, NPTS or CUTOFF have to be specified")
417 : END IF
418 :
419 : CALL pw_grid_setup_internal(cell_hmat, cell_h_inv, cell_deth, pw_grid, mp_comm, bounds_local=bounds_local, &
420 44594 : blocked=blocked, ref_grid=ref_grid, rs_dims=rs_dims, iounit=iounit)
421 :
422 44594 : CALL pw_grid_create_ghatmap(pw_grid)
423 :
424 44594 : CALL timestop(handle)
425 :
426 44594 : END SUBROUTINE pw_grid_create_extended
427 :
428 : ! **************************************************************************************************
429 : !> \brief sets up a combined index for gather and scatter
430 : !> \param pw_grid ...
431 : !> \author Gloess Andreas (xx-Dec-2012)
432 : ! **************************************************************************************************
433 44594 : SUBROUTINE pw_grid_create_ghatmap(pw_grid)
434 :
435 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
436 :
437 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_create_ghatmap'
438 :
439 : INTEGER :: gpt, handle, l, m, mn, n
440 :
441 44594 : CALL timeset(routineN, handle)
442 :
443 : ! some checks
444 44594 : CPASSERT(pw_grid%ref_count > 0)
445 :
446 : ! mapping of map_x( g_hat(i,j)) to g_hatmap
447 : ! the second index is for switching from gather(1) to scatter(2)
448 : ASSOCIATE (g_hat => pw_grid%g_hat, g_hatmap => pw_grid%g_hatmap, pmapl => pw_grid%mapl%pos, &
449 : pmapm => pw_grid%mapm%pos, pmapn => pw_grid%mapn%pos, nmapl => pw_grid%mapl%neg, &
450 : nmapm => pw_grid%mapm%neg, nmapn => pw_grid%mapn%neg, ngpts => SIZE(pw_grid%gsq), &
451 : npts => pw_grid%npts, yzq => pw_grid%para%yzq)
452 : ! initialize map array to minus one, to guarantee memory
453 : ! range checking errors (just to be sure)
454 1162961034 : g_hatmap(:, :) = -1
455 89188 : IF (pw_grid%para%mode /= PW_MODE_DISTRIBUTED) THEN
456 21212 : CPASSERT(PRODUCT(INT(npts, KIND=int_8)) <= INT(HUGE(0), KIND=int_8))
457 48663414 : DO gpt = 1, ngpts
458 48658111 : l = pmapl(g_hat(1, gpt))
459 48658111 : m = pmapm(g_hat(2, gpt))
460 48658111 : n = pmapn(g_hat(3, gpt))
461 : !ATTENTION: C-mapping [start-index=0] !!!!
462 : !ATTENTION: potential integer overflow !!!!
463 48663414 : g_hatmap(gpt, 1) = l + npts(1)*(m + npts(2)*n)
464 : END DO
465 5303 : IF (pw_grid%grid_span == HALFSPACE) THEN
466 12630608 : DO gpt = 1, ngpts
467 12629401 : l = nmapl(g_hat(1, gpt))
468 12629401 : m = nmapm(g_hat(2, gpt))
469 12629401 : n = nmapn(g_hat(3, gpt))
470 : !ATTENTION: C-mapping [start-index=0] !!!!
471 : !ATTENTION: potential integer overflow !!!!
472 12630608 : g_hatmap(gpt, 2) = l + npts(1)*(m + npts(2)*n)
473 : END DO
474 : END IF
475 : ELSE
476 1082905344 : DO gpt = 1, ngpts
477 1082866053 : l = pmapl(g_hat(1, gpt))
478 1082866053 : m = pmapm(g_hat(2, gpt)) + 1
479 1082866053 : n = pmapn(g_hat(3, gpt)) + 1
480 : !ATTENTION: C-mapping [start-index=0] !!!!
481 : !ATTENTION: potential integer overflow !!!!
482 1082866053 : mn = yzq(m, n) - 1
483 1082905344 : g_hatmap(gpt, 1) = l + npts(1)*mn
484 : END DO
485 39291 : IF (pw_grid%grid_span == HALFSPACE) THEN
486 18717074 : DO gpt = 1, ngpts
487 18714538 : l = nmapl(g_hat(1, gpt))
488 18714538 : m = nmapm(g_hat(2, gpt)) + 1
489 18714538 : n = nmapn(g_hat(3, gpt)) + 1
490 : !ATTENTION: C-mapping [start-index=0] !!!!
491 : !ATTENTION: potential integer overflow !!!!
492 18714538 : mn = yzq(m, n) - 1
493 18717074 : g_hatmap(gpt, 2) = l + npts(1)*mn
494 : END DO
495 : END IF
496 : END IF
497 : END ASSOCIATE
498 :
499 44594 : CALL timestop(handle)
500 :
501 44594 : END SUBROUTINE pw_grid_create_ghatmap
502 :
503 : ! **************************************************************************************************
504 : !> \brief sets up a pw_grid, needs valid bounds as input, it is up to you to
505 : !> make sure of it using pw_grid_bounds_from_n
506 : !> \param cell_hmat ...
507 : !> \param cell_h_inv ...
508 : !> \param cell_deth ...
509 : !> \param pw_grid ...
510 : !> \param mp_comm ...
511 : !> \param bounds_local ...
512 : !> \param blocked ...
513 : !> \param ref_grid ...
514 : !> \param rs_dims ...
515 : !> \param iounit ...
516 : !> \par History
517 : !> JGH (20-Dec-2000) : Adapted for parallel use
518 : !> JGH (28-Feb-2001) : New optional argument fft_usage
519 : !> JGH (21-Mar-2001) : Reference grid code
520 : !> JGH (21-Mar-2001) : New optional argument symm_usage
521 : !> JGH (22-Mar-2001) : Simplify group assignment (mp_comm_dup)
522 : !> JGH (21-May-2002) : Remove orthorhombic keyword (code is fast enough)
523 : !> JGH (19-Feb-2003) : Negative cutoff can be used for non-spheric grids
524 : !> Joost VandeVondele (Feb-2004) : optionally generate pw grids that are commensurate in rs
525 : !> JGH (18-Dec-2007) : Refactoring
526 : !> \author fawzi
527 : !> \note
528 : !> this is the function that should be used in the future
529 : ! **************************************************************************************************
530 44594 : SUBROUTINE pw_grid_setup_internal(cell_hmat, cell_h_inv, cell_deth, pw_grid, mp_comm, bounds_local, &
531 : blocked, ref_grid, rs_dims, iounit)
532 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: cell_hmat, cell_h_inv
533 : REAL(KIND=dp), INTENT(IN) :: cell_deth
534 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
535 :
536 : CLASS(mp_comm_type), INTENT(IN) :: mp_comm
537 : INTEGER, DIMENSION(2, 3), INTENT(IN), OPTIONAL :: bounds_local
538 : INTEGER, INTENT(in), OPTIONAL :: blocked
539 : TYPE(pw_grid_type), INTENT(in), OPTIONAL :: ref_grid
540 : INTEGER, DIMENSION(2), INTENT(in), OPTIONAL :: rs_dims
541 : INTEGER, INTENT(in), OPTIONAL :: iounit
542 :
543 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_setup_internal'
544 :
545 : INTEGER :: handle, n(3)
546 44594 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: yz_mask
547 : REAL(KIND=dp) :: ecut
548 :
549 : !------------------------------------------------------------------------------
550 :
551 44594 : CALL timeset(routineN, handle)
552 :
553 44594 : CPASSERT(pw_grid%ref_count > 0)
554 :
555 : ! set pointer to possible reference grid
556 44594 : IF (PRESENT(ref_grid)) THEN
557 27156 : pw_grid%reference = ref_grid%id_nr
558 : END IF
559 :
560 44594 : IF (pw_grid%spherical) THEN
561 3491 : ecut = pw_grid%cutoff
562 : ELSE
563 41103 : ecut = 1.e10_dp
564 : END IF
565 :
566 178376 : n(:) = pw_grid%npts(:)
567 :
568 : ! Find the number of grid points
569 : ! yz_mask counts the number of g-vectors orthogonal to the yz plane
570 : ! the indices in yz_mask are from -n/2 .. n/2 shifted by n/2 + 1
571 : ! these are not mapped indices !
572 178376 : ALLOCATE (yz_mask(n(2), n(3)))
573 44594 : CALL pw_grid_count(cell_h_inv, pw_grid, mp_comm, ecut, yz_mask)
574 :
575 : ! Check if reference grid is compatible
576 44594 : IF (PRESENT(ref_grid)) THEN
577 27156 : CPASSERT(pw_grid%para%mode == ref_grid%para%mode)
578 27156 : CPASSERT(pw_grid%grid_span == ref_grid%grid_span)
579 27156 : CPASSERT(pw_grid%spherical .EQV. ref_grid%spherical)
580 : END IF
581 :
582 : ! Distribute grid
583 : CALL pw_grid_distribute(pw_grid, mp_comm, yz_mask, bounds_local=bounds_local, ref_grid=ref_grid, blocked=blocked, &
584 44594 : rs_dims=rs_dims)
585 :
586 : ! Allocate the grid fields
587 : CALL pw_grid_allocate(pw_grid, pw_grid%ngpts_cut_local, &
588 44594 : pw_grid%bounds)
589 :
590 : ! Fill in the grid structure
591 44594 : CALL pw_grid_assign(cell_h_inv, pw_grid, ecut)
592 :
593 : ! Sort g vector wrt length (only local for each processor)
594 44594 : CALL pw_grid_sort(pw_grid, ref_grid)
595 :
596 44594 : CALL pw_grid_remap(pw_grid, yz_mask)
597 :
598 44594 : DEALLOCATE (yz_mask)
599 :
600 44594 : CALL cell2grid(cell_hmat, cell_h_inv, cell_deth, pw_grid)
601 : !
602 : ! Output: All the information of this grid type
603 : !
604 :
605 44594 : IF (PRESENT(iounit)) THEN
606 44568 : CALL pw_grid_print(pw_grid, iounit)
607 : END IF
608 :
609 44594 : CALL timestop(handle)
610 :
611 44594 : END SUBROUTINE pw_grid_setup_internal
612 :
613 : ! **************************************************************************************************
614 : !> \brief Helper routine used by pw_grid_setup_internal and pw_grid_change
615 : !> \param cell_hmat ...
616 : !> \param cell_h_inv ...
617 : !> \param cell_deth ...
618 : !> \param pw_grid ...
619 : !> \par History moved common code into new subroutine
620 : !> \author Ole Schuett
621 : ! **************************************************************************************************
622 58474 : SUBROUTINE cell2grid(cell_hmat, cell_h_inv, cell_deth, pw_grid)
623 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: cell_hmat, cell_h_inv
624 : REAL(KIND=dp), INTENT(IN) :: cell_deth
625 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
626 :
627 58474 : pw_grid%vol = ABS(cell_deth)
628 58474 : pw_grid%dvol = pw_grid%vol/REAL(pw_grid%ngpts, KIND=dp)
629 : pw_grid%dr(1) = SQRT(SUM(cell_hmat(:, 1)**2)) &
630 233896 : /REAL(pw_grid%npts(1), KIND=dp)
631 : pw_grid%dr(2) = SQRT(SUM(cell_hmat(:, 2)**2)) &
632 233896 : /REAL(pw_grid%npts(2), KIND=dp)
633 : pw_grid%dr(3) = SQRT(SUM(cell_hmat(:, 3)**2)) &
634 233896 : /REAL(pw_grid%npts(3), KIND=dp)
635 233896 : pw_grid%dh(:, 1) = cell_hmat(:, 1)/REAL(pw_grid%npts(1), KIND=dp)
636 233896 : pw_grid%dh(:, 2) = cell_hmat(:, 2)/REAL(pw_grid%npts(2), KIND=dp)
637 233896 : pw_grid%dh(:, 3) = cell_hmat(:, 3)/REAL(pw_grid%npts(3), KIND=dp)
638 233896 : pw_grid%dh_inv(1, :) = cell_h_inv(1, :)*REAL(pw_grid%npts(1), KIND=dp)
639 233896 : pw_grid%dh_inv(2, :) = cell_h_inv(2, :)*REAL(pw_grid%npts(2), KIND=dp)
640 233896 : pw_grid%dh_inv(3, :) = cell_h_inv(3, :)*REAL(pw_grid%npts(3), KIND=dp)
641 :
642 : IF ((cell_hmat(1, 2) == 0.0_dp) .AND. (cell_hmat(1, 3) == 0.0_dp) .AND. &
643 : (cell_hmat(2, 1) == 0.0_dp) .AND. (cell_hmat(2, 3) == 0.0_dp) .AND. &
644 58474 : (cell_hmat(3, 1) == 0.0_dp) .AND. (cell_hmat(3, 2) == 0.0_dp)) THEN
645 45904 : pw_grid%orthorhombic = .TRUE.
646 : ELSE
647 12570 : pw_grid%orthorhombic = .FALSE.
648 : END IF
649 58474 : END SUBROUTINE cell2grid
650 :
651 : ! **************************************************************************************************
652 : !> \brief Output of information on pw_grid
653 : !> \param pw_grid ...
654 : !> \param info ...
655 : !> \author JGH[18-05-2007] from earlier versions
656 : ! **************************************************************************************************
657 44568 : SUBROUTINE pw_grid_print(pw_grid, info)
658 :
659 : TYPE(pw_grid_type), INTENT(IN) :: pw_grid
660 : INTEGER, INTENT(IN) :: info
661 :
662 : INTEGER :: i
663 : INTEGER(KIND=int_8) :: n(3)
664 : REAL(KIND=dp) :: rv(3, 3)
665 :
666 : !------------------------------------------------------------------------------
667 : !
668 : ! Output: All the information of this grid type
669 : !
670 :
671 44568 : IF (pw_grid%para%mode == PW_MODE_LOCAL) THEN
672 5303 : IF (info > 0) THEN
673 : WRITE (info, '(/,A,T71,I10)') &
674 447 : " PW_GRID| Information for grid number ", pw_grid%id_nr
675 447 : IF (pw_grid%reference > 0) THEN
676 : WRITE (info, '(A,T71,I10)') &
677 114 : " PW_GRID| Number of the reference grid ", pw_grid%reference
678 : END IF
679 447 : WRITE (info, '(" PW_GRID| Cutoff [a.u.]",T71,f10.1)') pw_grid%cutoff
680 447 : IF (pw_grid%spherical) THEN
681 201 : WRITE (info, '(A,T78,A)') " PW_GRID| spherical cutoff: ", "YES"
682 201 : WRITE (info, '(A,T71,I10)') " PW_GRID| Grid points within cutoff", &
683 402 : pw_grid%ngpts_cut
684 : ELSE
685 246 : WRITE (info, '(A,T78,A)') " PW_GRID| spherical cutoff: ", " NO"
686 : END IF
687 1788 : DO i = 1, 3
688 1341 : WRITE (info, '(A,I3,T30,2I8,T62,A,T71,I10)') " PW_GRID| Bounds ", &
689 1341 : i, pw_grid%bounds(1, I), pw_grid%bounds(2, I), &
690 3129 : "Points:", pw_grid%npts(I)
691 : END DO
692 : WRITE (info, '(A,G12.4,T50,A,T67,F14.4)') &
693 447 : " PW_GRID| Volume element (a.u.^3)", &
694 894 : pw_grid%dvol, " Volume (a.u.^3) :", pw_grid%vol
695 447 : IF (pw_grid%grid_span == HALFSPACE) THEN
696 299 : WRITE (info, '(A,T72,A)') " PW_GRID| Grid span", "HALFSPACE"
697 : ELSE
698 148 : WRITE (info, '(A,T72,A)') " PW_GRID| Grid span", "FULLSPACE"
699 : END IF
700 : END IF
701 : ELSE
702 :
703 39265 : n(1) = pw_grid%ngpts_cut_local
704 39265 : n(2) = pw_grid%ngpts_local
705 39265 : CALL pw_grid%para%group%sum(n(1:2))
706 117795 : n(3) = SUM(pw_grid%para%nyzray)
707 157060 : rv(:, 1) = REAL(n, KIND=dp)/REAL(pw_grid%para%group%num_pe, KIND=dp)
708 39265 : n(1) = pw_grid%ngpts_cut_local
709 39265 : n(2) = pw_grid%ngpts_local
710 39265 : CALL pw_grid%para%group%max(n(1:2))
711 117795 : n(3) = MAXVAL(pw_grid%para%nyzray)
712 157060 : rv(:, 2) = REAL(n, KIND=dp)
713 39265 : n(1) = pw_grid%ngpts_cut_local
714 39265 : n(2) = pw_grid%ngpts_local
715 39265 : CALL pw_grid%para%group%min(n(1:2))
716 117795 : n(3) = MINVAL(pw_grid%para%nyzray)
717 157060 : rv(:, 3) = REAL(n, KIND=dp)
718 :
719 39265 : IF (info > 0) THEN
720 : WRITE (info, '(/,A,T71,I10)') &
721 8015 : " PW_GRID| Information for grid number ", pw_grid%id_nr
722 8015 : IF (pw_grid%reference > 0) THEN
723 : WRITE (info, '(A,T71,I10)') &
724 4799 : " PW_GRID| Number of the reference grid ", pw_grid%reference
725 : END IF
726 : WRITE (info, '(A,T60,I10,A)') &
727 8015 : " PW_GRID| Grid distributed over ", pw_grid%para%group%num_pe, &
728 16030 : " processors"
729 : WRITE (info, '(A,T71,2I5)') &
730 8015 : " PW_GRID| Real space group dimensions ", pw_grid%para%group%num_pe_cart
731 8015 : IF (pw_grid%para%blocked) THEN
732 6 : WRITE (info, '(A,T78,A)') " PW_GRID| the grid is blocked: ", "YES"
733 : ELSE
734 8009 : WRITE (info, '(A,T78,A)') " PW_GRID| the grid is blocked: ", " NO"
735 : END IF
736 8015 : WRITE (info, '(" PW_GRID| Cutoff [a.u.]",T71,f10.1)') pw_grid%cutoff
737 8015 : IF (pw_grid%spherical) THEN
738 408 : WRITE (info, '(A,T78,A)') " PW_GRID| spherical cutoff: ", "YES"
739 408 : WRITE (info, '(A,T71,I10)') " PW_GRID| Grid points within cutoff", &
740 816 : pw_grid%ngpts_cut
741 : ELSE
742 7607 : WRITE (info, '(A,T78,A)') " PW_GRID| spherical cutoff: ", " NO"
743 : END IF
744 32060 : DO i = 1, 3
745 24045 : WRITE (info, '(A,I3,T30,2I8,T62,A,T71,I10)') " PW_GRID| Bounds ", &
746 24045 : i, pw_grid%bounds(1, I), pw_grid%bounds(2, I), &
747 56105 : "Points:", pw_grid%npts(I)
748 : END DO
749 : WRITE (info, '(A,G12.4,T50,A,T67,F14.4)') &
750 8015 : " PW_GRID| Volume element (a.u.^3)", &
751 16030 : pw_grid%dvol, " Volume (a.u.^3) :", pw_grid%vol
752 8015 : IF (pw_grid%grid_span == HALFSPACE) THEN
753 416 : WRITE (info, '(A,T72,A)') " PW_GRID| Grid span", "HALFSPACE"
754 : ELSE
755 7599 : WRITE (info, '(A,T72,A)') " PW_GRID| Grid span", "FULLSPACE"
756 : END IF
757 8015 : WRITE (info, '(A,T48,A)') " PW_GRID| Distribution", &
758 16030 : " Average Max Min"
759 8015 : WRITE (info, '(A,T45,F12.1,2I12)') " PW_GRID| G-Vectors", &
760 16030 : rv(1, 1), NINT(rv(1, 2)), NINT(rv(1, 3))
761 8015 : WRITE (info, '(A,T45,F12.1,2I12)') " PW_GRID| G-Rays ", &
762 16030 : rv(3, 1), NINT(rv(3, 2)), NINT(rv(3, 3))
763 8015 : WRITE (info, '(A,T45,F12.1,2I12)') " PW_GRID| Real Space Points", &
764 16030 : rv(2, 1), NINT(rv(2, 2)), NINT(rv(2, 3))
765 : END IF ! group head
766 : END IF ! local
767 :
768 44568 : END SUBROUTINE pw_grid_print
769 :
770 : ! **************************************************************************************************
771 : !> \brief Distribute grids in real and Fourier Space to the processors in group
772 : !> \param pw_grid ...
773 : !> \param mp_comm ...
774 : !> \param yz_mask ...
775 : !> \param bounds_local ...
776 : !> \param ref_grid ...
777 : !> \param blocked ...
778 : !> \param rs_dims ...
779 : !> \par History
780 : !> JGH (01-Mar-2001) optional reference grid
781 : !> JGH (22-May-2002) bug fix for pre_tag and HALFSPACE grids
782 : !> JGH (09-Sep-2003) reduce scaling for distribution
783 : !> \author JGH (22-12-2000)
784 : ! **************************************************************************************************
785 44594 : SUBROUTINE pw_grid_distribute(pw_grid, mp_comm, yz_mask, bounds_local, ref_grid, blocked, rs_dims)
786 :
787 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
788 :
789 : CLASS(mp_comm_type), INTENT(IN) :: mp_comm
790 : INTEGER, DIMENSION(:, :), INTENT(INOUT) :: yz_mask
791 : INTEGER, DIMENSION(2, 3), INTENT(IN), OPTIONAL :: bounds_local
792 : TYPE(pw_grid_type), INTENT(IN), OPTIONAL :: ref_grid
793 : INTEGER, INTENT(IN), OPTIONAL :: blocked
794 : INTEGER, DIMENSION(2), INTENT(in), OPTIONAL :: rs_dims
795 :
796 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_distribute'
797 :
798 : INTEGER :: blocked_local, coor(2), gmax, handle, i, &
799 : i1, i2, ip, ipl, ipp, itmp, j, l, lby, &
800 : lbz, lo(2), m, n, np, ns, nx, ny, nz, &
801 : rsd(2)
802 44594 : INTEGER, ALLOCATABLE, DIMENSION(:) :: pemap
803 44594 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: yz_index
804 44594 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: axis_dist_all
805 : INTEGER, DIMENSION(2) :: my_rs_dims
806 : INTEGER, DIMENSION(2, 3) :: axis_dist
807 : LOGICAL :: blocking
808 :
809 : !------------------------------------------------------------------------------
810 :
811 44594 : CALL timeset(routineN, handle)
812 :
813 44594 : lby = pw_grid%bounds(1, 2)
814 44594 : lbz = pw_grid%bounds(1, 3)
815 :
816 178376 : pw_grid%ngpts = PRODUCT(INT(pw_grid%npts, KIND=int_8))
817 :
818 44594 : my_rs_dims = 0
819 44594 : IF (PRESENT(rs_dims)) THEN
820 43525 : my_rs_dims = rs_dims
821 : END IF
822 :
823 44594 : IF (PRESENT(blocked)) THEN
824 40879 : blocked_local = blocked
825 : ELSE
826 : blocked_local = do_pw_grid_blocked_free
827 : END IF
828 :
829 44594 : pw_grid%para%blocked = .FALSE.
830 :
831 44594 : IF (pw_grid%para%mode == PW_MODE_LOCAL) THEN
832 :
833 5303 : pw_grid%para%ray_distribution = .FALSE.
834 :
835 53030 : pw_grid%bounds_local = pw_grid%bounds
836 21212 : pw_grid%npts_local = pw_grid%npts
837 5303 : CPASSERT(pw_grid%ngpts_cut < HUGE(pw_grid%ngpts_cut_local))
838 5303 : pw_grid%ngpts_cut_local = INT(pw_grid%ngpts_cut)
839 5303 : CPASSERT(pw_grid%ngpts < HUGE(pw_grid%ngpts_local))
840 5303 : pw_grid%ngpts_local = INT(pw_grid%ngpts)
841 15909 : my_rs_dims = 1
842 5303 : CALL pw_grid%para%group%create(mp_comm, 2, my_rs_dims)
843 :
844 5303 : ALLOCATE (pw_grid%para%bo(2, 3, 0:0, 3))
845 21212 : DO i = 1, 3
846 63636 : pw_grid%para%bo(1, 1:3, 0, i) = 1
847 68939 : pw_grid%para%bo(2, 1:3, 0, i) = pw_grid%npts(1:3)
848 : END DO
849 :
850 : ELSE
851 :
852 : !..find the real space distribution
853 39291 : nx = pw_grid%npts(1)
854 39291 : ny = pw_grid%npts(2)
855 39291 : nz = pw_grid%npts(3)
856 39291 : np = mp_comm%num_pe
857 :
858 : ! The user can specify 2 strictly positive indices => specific layout
859 : ! 1 strictly positive index => the other is fixed by the number of CPUs
860 : ! 0 strictly positive indices => fully free distribution
861 : ! if fully free or the user request can not be fulfilled, we optimize heuristically ourselves by:
862 : ! 1) nx>np -> taking a plane distribution (/np,1/)
863 : ! 2) nx<np -> taking the most square distribution
864 : ! if blocking is free:
865 : ! 1) blocked=.FALSE. for plane distributions
866 : ! 2) blocked=.TRUE. for non-plane distributions
867 53131 : IF (ANY(my_rs_dims <= 0)) THEN
868 97089 : IF (ALL(my_rs_dims <= 0)) THEN
869 32347 : my_rs_dims = [0, 0]
870 : ELSE
871 24 : IF (my_rs_dims(1) > 0) THEN
872 0 : my_rs_dims(2) = np/my_rs_dims(1)
873 : ELSE
874 24 : my_rs_dims(1) = np/my_rs_dims(2)
875 : END IF
876 : END IF
877 : END IF
878 : ! reset if the distribution can not be fulfilled
879 117873 : IF (PRODUCT(my_rs_dims) /= np) my_rs_dims = [0, 0]
880 : ! reset if the distribution can not be dealt with [1,np]
881 39339 : IF (ALL(my_rs_dims == [1, np])) my_rs_dims = [0, 0]
882 :
883 : ! if [0,0] now, we can optimize it ourselves
884 104041 : IF (ALL(my_rs_dims == [0, 0])) THEN
885 : ! only small grids have a chance to be 2d distributed
886 32375 : IF (nx < np) THEN
887 : ! gives the most square looking distribution
888 0 : CALL mp_dims_create(np, my_rs_dims)
889 : ! we tend to like the first index being smaller than the second
890 0 : IF (my_rs_dims(1) > my_rs_dims(2)) THEN
891 0 : itmp = my_rs_dims(1)
892 0 : my_rs_dims(1) = my_rs_dims(2)
893 0 : my_rs_dims(2) = itmp
894 : END IF
895 : ! but should avoid having the first index 1 in all cases
896 0 : IF (my_rs_dims(1) == 1) THEN
897 0 : itmp = my_rs_dims(1)
898 0 : my_rs_dims(1) = my_rs_dims(2)
899 0 : my_rs_dims(2) = itmp
900 : END IF
901 : ELSE
902 97125 : my_rs_dims = [np, 1]
903 : END IF
904 : END IF
905 :
906 : ! now fix the blocking if we have a choice
907 : SELECT CASE (blocked_local)
908 : CASE (do_pw_grid_blocked_false)
909 28 : blocking = .FALSE.
910 : CASE (do_pw_grid_blocked_true)
911 28 : blocking = .TRUE.
912 : CASE (do_pw_grid_blocked_free)
913 40137 : IF (ALL(my_rs_dims == [np, 1])) THEN
914 : blocking = .FALSE.
915 : ELSE
916 0 : blocking = .TRUE.
917 : END IF
918 : CASE DEFAULT
919 39291 : CPABORT("Unknown PW grid blocking")
920 : END SELECT
921 :
922 : !..create group for real space distribution
923 39291 : CALL pw_grid%para%group%create(mp_comm, 2, my_rs_dims)
924 :
925 39291 : IF (PRESENT(bounds_local)) THEN
926 220 : pw_grid%bounds_local = bounds_local
927 : ELSE
928 : lo = get_limit(nx, pw_grid%para%group%num_pe_cart(1), &
929 39269 : pw_grid%para%group%mepos_cart(1))
930 117807 : pw_grid%bounds_local(:, 1) = lo + pw_grid%bounds(1, 1) - 1
931 : lo = get_limit(ny, pw_grid%para%group%num_pe_cart(2), &
932 39269 : pw_grid%para%group%mepos_cart(2))
933 117807 : pw_grid%bounds_local(:, 2) = lo + pw_grid%bounds(1, 2) - 1
934 117807 : pw_grid%bounds_local(:, 3) = pw_grid%bounds(:, 3)
935 : END IF
936 :
937 : pw_grid%npts_local(:) = pw_grid%bounds_local(2, :) &
938 157164 : - pw_grid%bounds_local(1, :) + 1
939 :
940 : !..the third distribution is needed for the second step in the FFT
941 157164 : ALLOCATE (pw_grid%para%bo(2, 3, 0:np - 1, 3))
942 117873 : rsd = pw_grid%para%group%num_pe_cart
943 :
944 39291 : IF (PRESENT(bounds_local)) THEN
945 : ! axis_dist tells what portion of 1 .. nx , 1 .. ny , 1 .. nz are in the current process
946 88 : DO i = 1, 3
947 220 : axis_dist(:, i) = bounds_local(:, i) - pw_grid%bounds(1, i) + 1
948 : END DO
949 66 : ALLOCATE (axis_dist_all(2, 3, np))
950 22 : CALL pw_grid%para%group%allgather(axis_dist, axis_dist_all)
951 66 : DO ip = 0, np - 1
952 44 : CALL pw_grid%para%group%coords(ip, coor)
953 : ! distribution xyZ
954 132 : pw_grid%para%bo(1:2, 1, ip, 1) = axis_dist_all(1:2, 1, ip + 1)
955 132 : pw_grid%para%bo(1:2, 2, ip, 1) = axis_dist_all(1:2, 2, ip + 1)
956 44 : pw_grid%para%bo(1, 3, ip, 1) = 1
957 44 : pw_grid%para%bo(2, 3, ip, 1) = nz
958 : ! distribution xYz
959 132 : pw_grid%para%bo(1:2, 1, ip, 2) = axis_dist_all(1:2, 1, ip + 1)
960 44 : pw_grid%para%bo(1, 2, ip, 2) = 1
961 44 : pw_grid%para%bo(2, 2, ip, 2) = ny
962 132 : pw_grid%para%bo(1:2, 3, ip, 2) = get_limit(nz, rsd(2), coor(2))
963 : ! distribution Xyz
964 44 : pw_grid%para%bo(1, 1, ip, 3) = 1
965 44 : pw_grid%para%bo(2, 1, ip, 3) = nx
966 132 : pw_grid%para%bo(1:2, 2, ip, 3) = get_limit(ny, rsd(1), coor(1))
967 154 : pw_grid%para%bo(1:2, 3, ip, 3) = get_limit(nz, rsd(2), coor(2))
968 : END DO
969 22 : DEALLOCATE (axis_dist_all)
970 : ELSE
971 117807 : DO ip = 0, np - 1
972 78538 : CALL pw_grid%para%group%coords(ip, coor)
973 : ! distribution xyZ
974 235614 : pw_grid%para%bo(1:2, 1, ip, 1) = get_limit(nx, rsd(1), coor(1))
975 235614 : pw_grid%para%bo(1:2, 2, ip, 1) = get_limit(ny, rsd(2), coor(2))
976 78538 : pw_grid%para%bo(1, 3, ip, 1) = 1
977 78538 : pw_grid%para%bo(2, 3, ip, 1) = nz
978 : ! distribution xYz
979 235614 : pw_grid%para%bo(1:2, 1, ip, 2) = get_limit(nx, rsd(1), coor(1))
980 78538 : pw_grid%para%bo(1, 2, ip, 2) = 1
981 78538 : pw_grid%para%bo(2, 2, ip, 2) = ny
982 235614 : pw_grid%para%bo(1:2, 3, ip, 2) = get_limit(nz, rsd(2), coor(2))
983 : ! distribution Xyz
984 78538 : pw_grid%para%bo(1, 1, ip, 3) = 1
985 78538 : pw_grid%para%bo(2, 1, ip, 3) = nx
986 235614 : pw_grid%para%bo(1:2, 2, ip, 3) = get_limit(ny, rsd(1), coor(1))
987 274883 : pw_grid%para%bo(1:2, 3, ip, 3) = get_limit(nz, rsd(2), coor(2))
988 : END DO
989 : END IF
990 : !..find the g space distribution
991 39291 : pw_grid%ngpts_cut_local = 0
992 :
993 117873 : ALLOCATE (pw_grid%para%nyzray(0:np - 1))
994 :
995 157164 : ALLOCATE (pw_grid%para%yzq(ny, nz))
996 :
997 : IF (pw_grid%spherical .OR. pw_grid%grid_span == HALFSPACE &
998 39291 : .OR. .NOT. blocking) THEN
999 :
1000 39279 : pw_grid%para%ray_distribution = .TRUE.
1001 :
1002 : ! the mirror tags below only handle lower-skewed or symmetric
1003 : ! frequency windows (the only ones the grid setup produces today)
1004 39279 : IF (pw_grid%grid_span == HALFSPACE) THEN
1005 2536 : CPASSERT(pw_grid%bounds(2, 2) <= -pw_grid%bounds(1, 2))
1006 2536 : CPASSERT(pw_grid%bounds(2, 3) <= -pw_grid%bounds(1, 3))
1007 : END IF
1008 :
1009 42025753 : pw_grid%para%yzq = -1
1010 39279 : IF (PRESENT(ref_grid)) THEN
1011 : ! tag all vectors from the reference grid
1012 24084 : CALL pre_tag(pw_grid, yz_mask, ref_grid)
1013 : END IF
1014 :
1015 : ! Round Robin distribution
1016 : ! Processors 0 .. NP-1, NP-1 .. 0 get the largest remaining batch
1017 : ! of g vectors in turn
1018 :
1019 39279 : i1 = SIZE(yz_mask, 1)
1020 39279 : i2 = SIZE(yz_mask, 2)
1021 117837 : ALLOCATE (yz_index(2, i1*i2))
1022 39279 : CALL order_mask(yz_mask, yz_index)
1023 41005840 : DO i = 1, i1*i2
1024 40966561 : lo(1) = yz_index(1, i)
1025 40966561 : lo(2) = yz_index(2, i)
1026 40966561 : IF (lo(1)*lo(2) == 0) CYCLE
1027 25534984 : gmax = yz_mask(lo(1), lo(2))
1028 25534984 : IF (gmax == 0) CYCLE
1029 25534984 : yz_mask(lo(1), lo(2)) = 0
1030 25534984 : ip = MOD(i - 1, 2*np)
1031 25534984 : IF (ip > np - 1) ip = 2*np - ip - 1
1032 25534984 : IF (ip == pw_grid%para%group%mepos) THEN
1033 12767492 : pw_grid%ngpts_cut_local = pw_grid%ngpts_cut_local + gmax
1034 : END IF
1035 25534984 : pw_grid%para%yzq(lo(1), lo(2)) = ip
1036 25574263 : IF (pw_grid%grid_span == HALFSPACE) THEN
1037 1166728 : m = -lo(1) - 2*lby + 2
1038 1166728 : n = -lo(2) - 2*lbz + 2
1039 1166728 : pw_grid%para%yzq(m, n) = ip
1040 1166728 : yz_mask(m, n) = 0
1041 : END IF
1042 : END DO
1043 :
1044 39279 : DEALLOCATE (yz_index)
1045 :
1046 : ! Count the total number of rays on each processor
1047 117837 : pw_grid%para%nyzray = 0
1048 1059192 : DO i = 1, nz
1049 42025753 : DO j = 1, ny
1050 40966561 : ip = pw_grid%para%yzq(j, i)
1051 40966561 : IF (ip >= 0) pw_grid%para%nyzray(ip) = &
1052 40812496 : pw_grid%para%nyzray(ip) + 1
1053 : END DO
1054 : END DO
1055 :
1056 : ! Allocate mapping array (y:z, nray, nproc)
1057 117837 : ns = MAXVAL(pw_grid%para%nyzray(0:np - 1))
1058 157116 : ALLOCATE (pw_grid%para%yzp(2, ns, 0:np - 1))
1059 :
1060 : ! Fill mapping array, recalculate nyzray for convenience
1061 117837 : pw_grid%para%nyzray = 0
1062 1059192 : DO i = 1, nz
1063 42025753 : DO j = 1, ny
1064 40966561 : ip = pw_grid%para%yzq(j, i)
1065 41986474 : IF (ip >= 0) THEN
1066 : pw_grid%para%nyzray(ip) = &
1067 39792583 : pw_grid%para%nyzray(ip) + 1
1068 39792583 : ns = pw_grid%para%nyzray(ip)
1069 39792583 : pw_grid%para%yzp(1, ns, ip) = j
1070 39792583 : pw_grid%para%yzp(2, ns, ip) = i
1071 39792583 : IF (ip == pw_grid%para%group%mepos) THEN
1072 19896304 : pw_grid%para%yzq(j, i) = ns
1073 : ELSE
1074 19896279 : pw_grid%para%yzq(j, i) = -1
1075 : END IF
1076 : ELSE
1077 1173978 : pw_grid%para%yzq(j, i) = -2
1078 : END IF
1079 : END DO
1080 : END DO
1081 :
1082 157116 : pw_grid%ngpts_local = PRODUCT(pw_grid%npts_local)
1083 :
1084 : ELSE
1085 : !
1086 : ! block distribution of g vectors, we do not have a spherical cutoff
1087 : !
1088 :
1089 12 : pw_grid%para%blocked = .TRUE.
1090 12 : pw_grid%para%ray_distribution = .FALSE.
1091 :
1092 36 : DO ip = 0, np - 1
1093 : m = pw_grid%para%bo(2, 2, ip, 3) - &
1094 24 : pw_grid%para%bo(1, 2, ip, 3) + 1
1095 : n = pw_grid%para%bo(2, 3, ip, 3) - &
1096 24 : pw_grid%para%bo(1, 3, ip, 3) + 1
1097 36 : pw_grid%para%nyzray(ip) = n*m
1098 : END DO
1099 :
1100 12 : ipl = pw_grid%para%group%mepos
1101 : l = pw_grid%para%bo(2, 1, ipl, 3) - &
1102 12 : pw_grid%para%bo(1, 1, ipl, 3) + 1
1103 : m = pw_grid%para%bo(2, 2, ipl, 3) - &
1104 12 : pw_grid%para%bo(1, 2, ipl, 3) + 1
1105 : n = pw_grid%para%bo(2, 3, ipl, 3) - &
1106 12 : pw_grid%para%bo(1, 3, ipl, 3) + 1
1107 12 : pw_grid%ngpts_cut_local = l*m*n
1108 12 : pw_grid%ngpts_local = pw_grid%ngpts_cut_local
1109 :
1110 4816 : pw_grid%para%yzq = 0
1111 : ny = pw_grid%para%bo(2, 2, ipl, 3) - &
1112 12 : pw_grid%para%bo(1, 2, ipl, 3) + 1
1113 254 : DO n = pw_grid%para%bo(1, 3, ipl, 3), &
1114 12 : pw_grid%para%bo(2, 3, ipl, 3)
1115 242 : i = n - pw_grid%para%bo(1, 3, ipl, 3)
1116 2523 : DO m = pw_grid%para%bo(1, 2, ipl, 3), &
1117 254 : pw_grid%para%bo(2, 2, ipl, 3)
1118 2281 : j = m - pw_grid%para%bo(1, 2, ipl, 3) + 1
1119 2523 : pw_grid%para%yzq(m, n) = j + i*ny
1120 : END DO
1121 : END DO
1122 :
1123 : ! Allocate mapping array (y:z, nray, nproc)
1124 36 : ns = MAXVAL(pw_grid%para%nyzray(0:np - 1))
1125 48 : ALLOCATE (pw_grid%para%yzp(2, ns, 0:np - 1))
1126 14112 : pw_grid%para%yzp = 0
1127 :
1128 24 : ALLOCATE (pemap(0:np - 1))
1129 12 : pemap = 0
1130 12 : pemap(pw_grid%para%group%mepos) = pw_grid%para%group%mepos
1131 12 : CALL pw_grid%para%group%sum(pemap)
1132 :
1133 36 : DO ip = 0, np - 1
1134 24 : ipp = pemap(ip)
1135 24 : ns = 0
1136 508 : DO n = pw_grid%para%bo(1, 3, ipp, 3), &
1137 24 : pw_grid%para%bo(2, 3, ipp, 3)
1138 484 : i = n - pw_grid%bounds(1, 3) + 1
1139 5046 : DO m = pw_grid%para%bo(1, 2, ipp, 3), &
1140 508 : pw_grid%para%bo(2, 2, ipp, 3)
1141 4562 : j = m - pw_grid%bounds(1, 2) + 1
1142 4562 : ns = ns + 1
1143 4562 : pw_grid%para%yzp(1, ns, ip) = j
1144 5046 : pw_grid%para%yzp(2, ns, ip) = i
1145 : END DO
1146 : END DO
1147 36 : CPASSERT(ns == pw_grid%para%nyzray(ip))
1148 : END DO
1149 :
1150 12 : DEALLOCATE (pemap)
1151 :
1152 : END IF
1153 :
1154 : END IF
1155 :
1156 : ! pos_of_x(i) tells on which cpu pw%array(i,:,:) is located
1157 : ! should be computable in principle, without the need for communication
1158 44594 : IF (pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
1159 117873 : ALLOCATE (pw_grid%para%pos_of_x(pw_grid%bounds(1, 1):pw_grid%bounds(2, 1)))
1160 1027678 : pw_grid%para%pos_of_x = 0
1161 533487 : pw_grid%para%pos_of_x(pw_grid%bounds_local(1, 1):pw_grid%bounds_local(2, 1)) = pw_grid%para%group%mepos
1162 39291 : CALL pw_grid%para%group%sum(pw_grid%para%pos_of_x)
1163 : ELSE
1164 : ! this should not be needed
1165 15909 : ALLOCATE (pw_grid%para%pos_of_x(pw_grid%bounds(1, 1):pw_grid%bounds(2, 1)))
1166 91704 : pw_grid%para%pos_of_x = 0
1167 : END IF
1168 :
1169 44594 : CALL timestop(handle)
1170 :
1171 44594 : END SUBROUTINE pw_grid_distribute
1172 :
1173 : ! **************************************************************************************************
1174 : !> \brief ...
1175 : !> \param pw_grid ...
1176 : !> \param yz_mask ...
1177 : !> \param ref_grid ...
1178 : !> \par History
1179 : !> - Fix mapping bug for pw_grid eqv to ref_grid (21.11.2019, MK)
1180 : ! **************************************************************************************************
1181 24084 : SUBROUTINE pre_tag(pw_grid, yz_mask, ref_grid)
1182 :
1183 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
1184 : INTEGER, DIMENSION(:, :), INTENT(INOUT) :: yz_mask
1185 : TYPE(pw_grid_type), INTENT(IN) :: ref_grid
1186 :
1187 : INTEGER :: gmax, ig, ip, lby, lbz, my, mz, ny, nz, &
1188 : uby, ubz, y, yp, z, zp
1189 :
1190 24084 : ny = ref_grid%npts(2)
1191 24084 : nz = ref_grid%npts(3)
1192 24084 : lby = pw_grid%bounds(1, 2)
1193 24084 : lbz = pw_grid%bounds(1, 3)
1194 24084 : uby = pw_grid%bounds(2, 2)
1195 24084 : ubz = pw_grid%bounds(2, 3)
1196 24084 : my = SIZE(yz_mask, 1)
1197 24084 : mz = SIZE(yz_mask, 2)
1198 :
1199 : ! loop over all processors and all g vectors yz lines on this processor
1200 72252 : DO ip = 0, ref_grid%para%group%num_pe - 1
1201 64604246 : DO ig = 1, ref_grid%para%nyzray(ip)
1202 : ! go from mapped coordinates to original coordinates
1203 : ! 1, 2, ..., n-1, n -> 0, 1, ..., (n/2)-1, -(n/2), -(n/2)+1, ..., -2, -1
1204 64531994 : y = ref_grid%para%yzp(1, ig, ip) - 1
1205 64531994 : IF (y >= ny/2) y = y - ny
1206 64531994 : z = ref_grid%para%yzp(2, ig, ip) - 1
1207 64531994 : IF (z >= nz/2) z = z - nz
1208 : ! check if this is inside the realm of the new grid
1209 64531994 : IF (y < lby .OR. y > uby .OR. z < lbz .OR. z > ubz) CYCLE
1210 : ! go to shifted coordinates
1211 13093407 : y = y - lby + 1
1212 13093407 : z = z - lbz + 1
1213 : ! this tag is outside the cutoff range of the new grid
1214 13093407 : IF (pw_grid%grid_span == HALFSPACE) THEN
1215 0 : yp = -y - 2*lby + 2
1216 0 : zp = -z - 2*lbz + 2
1217 : ! if the reference grid is larger than the mirror point may be
1218 : ! outside the new grid even if the original point is inside
1219 0 : IF (yp < 1 .OR. yp > my .OR. zp < 1 .OR. zp > mz) CYCLE
1220 0 : gmax = MAX(yz_mask(y, z), yz_mask(yp, zp))
1221 0 : IF (gmax == 0) CYCLE
1222 0 : yz_mask(y, z) = 0
1223 0 : yz_mask(yp, zp) = 0
1224 0 : pw_grid%para%yzq(y, z) = ip
1225 0 : pw_grid%para%yzq(yp, zp) = ip
1226 : ELSE
1227 13093407 : gmax = yz_mask(y, z)
1228 13093407 : IF (gmax == 0) CYCLE
1229 13093407 : yz_mask(y, z) = 0
1230 13093407 : pw_grid%para%yzq(y, z) = ip
1231 : END IF
1232 13141575 : IF (ip == pw_grid%para%group%mepos) THEN
1233 6546716 : pw_grid%ngpts_cut_local = pw_grid%ngpts_cut_local + gmax
1234 : END IF
1235 : END DO
1236 : END DO
1237 :
1238 24084 : END SUBROUTINE pre_tag
1239 :
1240 : ! **************************************************************************************************
1241 : !> \brief ...
1242 : !> \param yz_mask ...
1243 : !> \param yz_index ...
1244 : ! **************************************************************************************************
1245 39279 : PURE SUBROUTINE order_mask(yz_mask, yz_index)
1246 :
1247 : INTEGER, DIMENSION(:, :), INTENT(IN) :: yz_mask
1248 : INTEGER, DIMENSION(:, :), INTENT(OUT) :: yz_index
1249 :
1250 : INTEGER :: i1, i2, ic, icount, ii, im, jc, jj
1251 :
1252 : !NB load balance
1253 : !------------------------------------------------------------------------------
1254 : !NB spiral out from origin, so that even if overall grid is full and
1255 : !NB block distributed, spherical cutoff still leads to good load
1256 : !NB balance in cp_ddapc_apply_CD
1257 :
1258 39279 : i1 = SIZE(yz_mask, 1)
1259 39279 : i2 = SIZE(yz_mask, 2)
1260 122938962 : yz_index = 0
1261 :
1262 39279 : icount = 1
1263 39279 : ic = i1/2
1264 39279 : jc = i2/2
1265 39279 : ii = ic
1266 39279 : jj = jc
1267 39279 : IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1268 39279 : IF (yz_mask(ii, jj) /= 0) THEN
1269 15195 : yz_index(1, icount) = ii
1270 15195 : yz_index(2, icount) = jj
1271 15195 : icount = icount + 1
1272 : END IF
1273 : END IF
1274 584860 : DO im = 1, MAX(ic + 1, jc + 1)
1275 545581 : ii = ic - im
1276 14091920 : DO jj = jc - im, jc + im
1277 14091920 : IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1278 9728019 : IF (yz_mask(ii, jj) /= 0) THEN
1279 6127407 : yz_index(1, icount) = ii
1280 6127407 : yz_index(2, icount) = jj
1281 6127407 : icount = icount + 1
1282 : END IF
1283 : END IF
1284 : END DO
1285 545581 : ii = ic + im
1286 14091920 : DO jj = jc - im, jc + im
1287 14091920 : IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1288 10984934 : IF (yz_mask(ii, jj) /= 0) THEN
1289 6682807 : yz_index(1, icount) = ii
1290 6682807 : yz_index(2, icount) = jj
1291 6682807 : icount = icount + 1
1292 : END IF
1293 : END IF
1294 : END DO
1295 545581 : jj = jc - im
1296 13000758 : DO ii = ic - im + 1, ic + im - 1
1297 13000758 : IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1298 9504917 : IF (yz_mask(ii, jj) /= 0) THEN
1299 6349681 : yz_index(1, icount) = ii
1300 6349681 : yz_index(2, icount) = jj
1301 6349681 : icount = icount + 1
1302 : END IF
1303 : END IF
1304 : END DO
1305 13000758 : jj = jc + im
1306 13040037 : DO ii = ic - im + 1, ic + im - 1
1307 13000758 : IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1308 10709412 : IF (yz_mask(ii, jj) /= 0) THEN
1309 6359894 : yz_index(1, icount) = ii
1310 6359894 : yz_index(2, icount) = jj
1311 6359894 : icount = icount + 1
1312 : END IF
1313 : END IF
1314 : END DO
1315 : END DO
1316 :
1317 39279 : END SUBROUTINE order_mask
1318 : ! **************************************************************************************************
1319 : !> \brief compute the length of g vectors
1320 : !> \param h_inv ...
1321 : !> \param length_x ...
1322 : !> \param length_y ...
1323 : !> \param length_z ...
1324 : !> \param length ...
1325 : !> \param l ...
1326 : !> \param m ...
1327 : !> \param n ...
1328 : ! **************************************************************************************************
1329 2324628995 : PURE SUBROUTINE pw_vec_length(h_inv, length_x, length_y, length_z, length, l, m, n)
1330 :
1331 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: h_inv
1332 : REAL(KIND=dp), INTENT(OUT) :: length_x, length_y, length_z, length
1333 : INTEGER, INTENT(IN) :: l, m, n
1334 :
1335 : length_x &
1336 : = REAL(l, dp)*h_inv(1, 1) &
1337 : + REAL(m, dp)*h_inv(2, 1) &
1338 2324628995 : + REAL(n, dp)*h_inv(3, 1)
1339 : length_y &
1340 : = REAL(l, dp)*h_inv(1, 2) &
1341 : + REAL(m, dp)*h_inv(2, 2) &
1342 2324628995 : + REAL(n, dp)*h_inv(3, 2)
1343 : length_z &
1344 : = REAL(l, dp)*h_inv(1, 3) &
1345 : + REAL(m, dp)*h_inv(2, 3) &
1346 2324628995 : + REAL(n, dp)*h_inv(3, 3)
1347 :
1348 : ! enforce strict zero-ness in this case (compiler optimization)
1349 2324628995 : IF (l == 0 .AND. m == 0 .AND. n == 0) THEN
1350 49877 : length_x = 0
1351 49877 : length_y = 0
1352 49877 : length_z = 0
1353 : END IF
1354 :
1355 2324628995 : length_x = length_x*twopi
1356 2324628995 : length_y = length_y*twopi
1357 2324628995 : length_z = length_z*twopi
1358 :
1359 2324628995 : length = length_x**2 + length_y**2 + length_z**2
1360 :
1361 2324628995 : END SUBROUTINE pw_vec_length
1362 :
1363 : ! **************************************************************************************************
1364 : !> \brief True for rays that must stay out of the stored half-space: on an
1365 : !> asymmetric (even) frequency window the -g mirror of a ray on the
1366 : !> lower extreme plane falls outside the grid, while the ray
1367 : !> distribution and the negation maps both require the stored
1368 : !> half-space to be closed under negation (the z = 0 plane cut in
1369 : !> pw_grid_count and the tally check in pw_grid_remap serve the same
1370 : !> contract). Symmetric (odd) windows are never affected.
1371 : !> \param grid_span ...
1372 : !> \param m ray frequency in direction 2
1373 : !> \param n ray frequency in direction 3
1374 : !> \param bounds global frequency bounds of the grid
1375 : !> \return whether the -g mirror of the ray falls outside the grid
1376 : ! **************************************************************************************************
1377 1251745019 : PURE FUNCTION mirror_outside_bounds(grid_span, m, n, bounds) RESULT(outside)
1378 :
1379 : INTEGER, INTENT(IN) :: grid_span, m, n
1380 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bounds
1381 : LOGICAL :: outside
1382 :
1383 : outside = grid_span == HALFSPACE .AND. &
1384 : ((m == bounds(1, 2) .AND. bounds(2, 2) < -bounds(1, 2)) .OR. &
1385 1251745019 : (n == bounds(1, 3) .AND. bounds(2, 3) < -bounds(1, 3)))
1386 :
1387 1251745019 : END FUNCTION mirror_outside_bounds
1388 :
1389 : ! **************************************************************************************************
1390 : !> \brief Count total number of g vectors
1391 : !> \param h_inv ...
1392 : !> \param pw_grid ...
1393 : !> \param mp_comm ...
1394 : !> \param cutoff ...
1395 : !> \param yz_mask ...
1396 : !> \par History
1397 : !> JGH (22-12-2000) : Adapted for parallel use
1398 : !> \author apsi
1399 : !> Christopher Mundy
1400 : ! **************************************************************************************************
1401 44594 : SUBROUTINE pw_grid_count(h_inv, pw_grid, mp_comm, cutoff, yz_mask)
1402 :
1403 : REAL(KIND=dp), DIMENSION(3, 3) :: h_inv
1404 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
1405 :
1406 : CLASS(mp_comm_type), INTENT(IN) :: mp_comm
1407 : REAL(KIND=dp), INTENT(IN) :: cutoff
1408 : INTEGER, DIMENSION(:, :), INTENT(OUT) :: yz_mask
1409 :
1410 : INTEGER :: l, m, mm, n, n_upperlimit, nlim(2), nn
1411 : INTEGER(KIND=int_8) :: gpt
1412 : REAL(KIND=dp) :: length, length_x, length_y, length_z
1413 :
1414 : ASSOCIATE (bounds => pw_grid%bounds)
1415 :
1416 44594 : IF (pw_grid%grid_span == HALFSPACE) THEN
1417 : n_upperlimit = 0
1418 40851 : ELSE IF (pw_grid%grid_span == FULLSPACE) THEN
1419 40851 : n_upperlimit = bounds(2, 3)
1420 : ELSE
1421 0 : CPABORT("No type set for the grid")
1422 : END IF
1423 :
1424 : ! finds valid g-points within grid
1425 44594 : gpt = 0
1426 44594 : IF (pw_grid%para%mode == PW_MODE_LOCAL) THEN
1427 5303 : nlim(1) = bounds(1, 3)
1428 5303 : nlim(2) = n_upperlimit
1429 39291 : ELSE IF (pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
1430 39291 : n = n_upperlimit - bounds(1, 3) + 1
1431 39291 : nlim = get_limit(n, mp_comm%num_pe, mp_comm%mepos)
1432 117873 : nlim = nlim + bounds(1, 3) - 1
1433 : ELSE
1434 0 : CPABORT("para % mode not specified")
1435 : END IF
1436 :
1437 44529362 : yz_mask = 0
1438 646931 : DO n = nlim(1), nlim(2)
1439 557743 : nn = n - bounds(1, 3) + 1
1440 21857187 : DO m = bounds(1, 2), bounds(2, 2)
1441 21254850 : mm = m - bounds(1, 2) + 1
1442 1193597221 : DO l = bounds(1, 1), bounds(2, 1)
1443 1171784628 : IF (pw_grid%grid_span == HALFSPACE .AND. n == 0) THEN
1444 3003142 : IF ((m == 0 .AND. l > 0) .OR. (m > 0)) CYCLE
1445 : END IF
1446 : ! keep the stored half-space closed under negation: the ray
1447 : ! distribution and the negation maps rely on it (see also the
1448 : ! z = 0 plane cut above and mirror_outside_bounds)
1449 1170121149 : IF (mirror_outside_bounds(pw_grid%grid_span, m, n, bounds)) CYCLE
1450 :
1451 1168715295 : CALL pw_vec_length(h_inv, length_x, length_y, length_z, length, l, m, n)
1452 :
1453 1189970145 : IF (0.5_dp*length <= cutoff) THEN
1454 1131524529 : gpt = gpt + 1
1455 1131524529 : yz_mask(mm, nn) = yz_mask(mm, nn) + 1
1456 : END IF
1457 :
1458 : END DO
1459 : END DO
1460 : END DO
1461 : END ASSOCIATE
1462 :
1463 : ! number of g-vectors for grid
1464 44594 : IF (pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
1465 39291 : CALL mp_comm%sum(gpt)
1466 84021847 : CALL mp_comm%sum(yz_mask)
1467 : END IF
1468 44594 : pw_grid%ngpts_cut = gpt
1469 :
1470 44594 : END SUBROUTINE pw_grid_count
1471 :
1472 : ! **************************************************************************************************
1473 : !> \brief Setup maps from 1d to 3d space
1474 : !> \param h_inv ...
1475 : !> \param pw_grid ...
1476 : !> \param cutoff ...
1477 : !> \par History
1478 : !> JGH (29-12-2000) : Adapted for parallel use
1479 : !> \author apsi
1480 : !> Christopher Mundy
1481 : ! **************************************************************************************************
1482 44594 : SUBROUTINE pw_grid_assign(h_inv, pw_grid, cutoff)
1483 :
1484 : REAL(KIND=dp), DIMENSION(3, 3) :: h_inv
1485 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
1486 : REAL(KIND=dp), INTENT(IN) :: cutoff
1487 :
1488 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_assign'
1489 :
1490 : INTEGER :: gpt, handle, i, ip, l, lby, lbz, ll, m, &
1491 : mm, n, n_upperlimit, nn
1492 : INTEGER(KIND=int_8) :: gpt_global
1493 : INTEGER, DIMENSION(2, 3) :: bol, bounds
1494 : REAL(KIND=dp) :: length, length_x, length_y, length_z
1495 :
1496 44594 : CALL timeset(routineN, handle)
1497 :
1498 445940 : bounds = pw_grid%bounds
1499 44594 : lby = pw_grid%bounds(1, 2)
1500 44594 : lbz = pw_grid%bounds(1, 3)
1501 :
1502 44594 : IF (pw_grid%grid_span == HALFSPACE) THEN
1503 : n_upperlimit = 0
1504 40851 : ELSE IF (pw_grid%grid_span == FULLSPACE) THEN
1505 40851 : n_upperlimit = bounds(2, 3)
1506 : ELSE
1507 0 : CPABORT("No type set for the grid")
1508 : END IF
1509 :
1510 : ! finds valid g-points within grid
1511 44594 : IF (pw_grid%para%mode == PW_MODE_LOCAL) THEN
1512 5303 : gpt = 0
1513 73233 : DO n = bounds(1, 3), n_upperlimit
1514 1686879 : DO m = bounds(1, 2), bounds(2, 2)
1515 64782208 : DO l = bounds(1, 1), bounds(2, 1)
1516 63100632 : IF (pw_grid%grid_span == HALFSPACE .AND. n == 0) THEN
1517 1270369 : IF ((m == 0 .AND. l > 0) .OR. (m > 0)) CYCLE
1518 : END IF
1519 : ! keep the stored half-space closed under negation (see pw_grid_count)
1520 62291818 : IF (mirror_outside_bounds(pw_grid%grid_span, m, n, bounds)) CYCLE
1521 :
1522 62167234 : CALL pw_vec_length(h_inv, length_x, length_y, length_z, length, l, m, n)
1523 :
1524 63780880 : IF (0.5_dp*length <= cutoff) THEN
1525 48658111 : gpt = gpt + 1
1526 48658111 : pw_grid%g(1, gpt) = length_x
1527 48658111 : pw_grid%g(2, gpt) = length_y
1528 48658111 : pw_grid%g(3, gpt) = length_z
1529 48658111 : pw_grid%gsq(gpt) = length
1530 48658111 : pw_grid%g_hat(1, gpt) = l
1531 48658111 : pw_grid%g_hat(2, gpt) = m
1532 48658111 : pw_grid%g_hat(3, gpt) = n
1533 : END IF
1534 :
1535 : END DO
1536 : END DO
1537 : END DO
1538 :
1539 : ELSE
1540 :
1541 39291 : IF (pw_grid%para%ray_distribution) THEN
1542 :
1543 39279 : gpt = 0
1544 39279 : ip = pw_grid%para%group%mepos
1545 19935583 : DO i = 1, pw_grid%para%nyzray(ip)
1546 19896304 : n = pw_grid%para%yzp(2, i, ip) + lbz - 1
1547 19896304 : m = pw_grid%para%yzp(1, i, ip) + lby - 1
1548 19896304 : IF (n > n_upperlimit) CYCLE
1549 : ! keep the stored half-space closed under negation (see pw_grid_count)
1550 19332052 : IF (mirror_outside_bounds(pw_grid%grid_span, m, n, bounds)) CYCLE
1551 1113850363 : DO l = bounds(1, 1), bounds(2, 1)
1552 1094479032 : IF (pw_grid%grid_span == HALFSPACE .AND. n == 0) THEN
1553 1552935 : IF ((m == 0 .AND. l > 0) .OR. (m > 0)) CYCLE
1554 : END IF
1555 :
1556 1093703545 : CALL pw_vec_length(h_inv, length_x, length_y, length_z, length, l, m, n)
1557 :
1558 1113599849 : IF (0.5_dp*length <= cutoff) THEN
1559 1082823132 : gpt = gpt + 1
1560 1082823132 : pw_grid%g(1, gpt) = length_x
1561 1082823132 : pw_grid%g(2, gpt) = length_y
1562 1082823132 : pw_grid%g(3, gpt) = length_z
1563 1082823132 : pw_grid%gsq(gpt) = length
1564 1082823132 : pw_grid%g_hat(1, gpt) = l
1565 1082823132 : pw_grid%g_hat(2, gpt) = m
1566 1082823132 : pw_grid%g_hat(3, gpt) = n
1567 : END IF
1568 :
1569 : END DO
1570 : END DO
1571 :
1572 : ELSE
1573 :
1574 120 : bol = pw_grid%para%bo(:, :, pw_grid%para%group%mepos, 3)
1575 12 : gpt = 0
1576 254 : DO n = bounds(1, 3), bounds(2, 3)
1577 242 : IF (n < 0) THEN
1578 120 : nn = n + pw_grid%npts(3) + 1
1579 : ELSE
1580 122 : nn = n + 1
1581 : END IF
1582 242 : IF (nn < bol(1, 3) .OR. nn > bol(2, 3)) CYCLE
1583 4816 : DO m = bounds(1, 2), bounds(2, 2)
1584 4562 : IF (m < 0) THEN
1585 2216 : mm = m + pw_grid%npts(2) + 1
1586 : ELSE
1587 2346 : mm = m + 1
1588 : END IF
1589 4562 : IF (mm < bol(1, 2) .OR. mm > bol(2, 2)) CYCLE
1590 45444 : DO l = bounds(1, 1), bounds(2, 1)
1591 42921 : IF (l < 0) THEN
1592 21148 : ll = l + pw_grid%npts(1) + 1
1593 : ELSE
1594 21773 : ll = l + 1
1595 : END IF
1596 42921 : IF (ll < bol(1, 1) .OR. ll > bol(2, 1)) CYCLE
1597 :
1598 42921 : CALL pw_vec_length(h_inv, length_x, length_y, length_z, length, l, m, n)
1599 :
1600 42921 : gpt = gpt + 1
1601 42921 : pw_grid%g(1, gpt) = length_x
1602 42921 : pw_grid%g(2, gpt) = length_y
1603 42921 : pw_grid%g(3, gpt) = length_z
1604 42921 : pw_grid%gsq(gpt) = length
1605 42921 : pw_grid%g_hat(1, gpt) = l
1606 42921 : pw_grid%g_hat(2, gpt) = m
1607 47483 : pw_grid%g_hat(3, gpt) = n
1608 :
1609 : END DO
1610 : END DO
1611 : END DO
1612 :
1613 : END IF
1614 :
1615 : END IF
1616 :
1617 : ! Check the number of g-vectors for grid
1618 44594 : CPASSERT(pw_grid%ngpts_cut_local == gpt)
1619 44594 : IF (pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
1620 39291 : gpt_global = gpt
1621 39291 : CALL pw_grid%para%group%sum(gpt_global)
1622 39291 : CPASSERT(pw_grid%ngpts_cut == gpt_global)
1623 : END IF
1624 :
1625 44594 : pw_grid%have_g0 = .FALSE.
1626 44594 : pw_grid%first_gne0 = 1
1627 851034970 : DO gpt = 1, pw_grid%ngpts_cut_local
1628 867125292 : IF (ALL(pw_grid%g_hat(:, gpt) == 0)) THEN
1629 24936 : pw_grid%have_g0 = .TRUE.
1630 24936 : pw_grid%first_gne0 = 2
1631 24936 : EXIT
1632 : END IF
1633 : END DO
1634 :
1635 : CALL pw_grid_set_maps(pw_grid%grid_span, pw_grid%g_hat, &
1636 44594 : pw_grid%mapl, pw_grid%mapm, pw_grid%mapn, pw_grid%npts)
1637 :
1638 44594 : CALL timestop(handle)
1639 :
1640 44594 : END SUBROUTINE pw_grid_assign
1641 :
1642 : ! **************************************************************************************************
1643 : !> \brief Setup maps from 1d to 3d space
1644 : !> \param grid_span ...
1645 : !> \param g_hat ...
1646 : !> \param mapl ...
1647 : !> \param mapm ...
1648 : !> \param mapn ...
1649 : !> \param npts ...
1650 : !> \par History
1651 : !> JGH (21-12-2000) : Size of g_hat locally determined
1652 : !> \author apsi
1653 : !> Christopher Mundy
1654 : !> \note
1655 : !> Maps are to full 3D space (not distributed)
1656 : ! **************************************************************************************************
1657 44594 : SUBROUTINE pw_grid_set_maps(grid_span, g_hat, mapl, mapm, mapn, npts)
1658 :
1659 : INTEGER, INTENT(IN) :: grid_span
1660 : INTEGER, DIMENSION(:, :), INTENT(IN) :: g_hat
1661 : TYPE(map_pn), INTENT(INOUT) :: mapl, mapm, mapn
1662 : INTEGER, DIMENSION(:), INTENT(IN) :: npts
1663 :
1664 : INTEGER :: gpt, l, m, n, ng
1665 :
1666 44594 : ng = SIZE(g_hat, 2)
1667 :
1668 1131568758 : DO gpt = 1, ng
1669 1131524164 : l = g_hat(1, gpt)
1670 1131524164 : m = g_hat(2, gpt)
1671 1131524164 : n = g_hat(3, gpt)
1672 1131524164 : IF (l < 0) THEN
1673 562665249 : mapl%pos(l) = l + npts(1)
1674 : ELSE
1675 568858915 : mapl%pos(l) = l
1676 : END IF
1677 1131524164 : IF (m < 0) THEN
1678 563150151 : mapm%pos(m) = m + npts(2)
1679 : ELSE
1680 568374013 : mapm%pos(m) = m
1681 : END IF
1682 1131524164 : IF (n < 0) THEN
1683 577983695 : mapn%pos(n) = n + npts(3)
1684 : ELSE
1685 553540469 : mapn%pos(n) = n
1686 : END IF
1687 :
1688 : ! Generating the maps to the full 3-d space
1689 :
1690 1131568758 : IF (grid_span == HALFSPACE) THEN
1691 :
1692 31343939 : IF (l <= 0) THEN
1693 16153597 : mapl%neg(l) = -l
1694 : ELSE
1695 15190342 : mapl%neg(l) = npts(1) - l
1696 : END IF
1697 31343939 : IF (m <= 0) THEN
1698 16575058 : mapm%neg(m) = -m
1699 : ELSE
1700 14768881 : mapm%neg(m) = npts(2) - m
1701 : END IF
1702 31343939 : IF (n <= 0) THEN
1703 31343939 : mapn%neg(n) = -n
1704 : ELSE
1705 0 : mapn%neg(n) = npts(3) - n
1706 : END IF
1707 :
1708 : END IF
1709 :
1710 : END DO
1711 :
1712 44594 : END SUBROUTINE pw_grid_set_maps
1713 :
1714 : ! **************************************************************************************************
1715 : !> \brief Allocate all (Pointer) Arrays in pw_grid
1716 : !> \param pw_grid ...
1717 : !> \param ng ...
1718 : !> \param bounds ...
1719 : !> \par History
1720 : !> JGH (20-12-2000) : Added status variable
1721 : !> Bounds of arrays now from calling routine, this
1722 : !> makes it independent from parallel setup
1723 : !> \author apsi
1724 : !> Christopher Mundy
1725 : ! **************************************************************************************************
1726 44594 : SUBROUTINE pw_grid_allocate(pw_grid, ng, bounds)
1727 :
1728 : ! Argument
1729 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
1730 : INTEGER, INTENT(IN) :: ng
1731 : INTEGER, DIMENSION(:, :), INTENT(IN) :: bounds
1732 :
1733 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_allocate'
1734 :
1735 : INTEGER :: nmaps
1736 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
1737 : INTEGER(KIND=C_SIZE_T) :: length
1738 : TYPE(C_PTR) :: cptr_g_hatmap
1739 : INTEGER :: stat
1740 : #endif
1741 :
1742 : INTEGER :: handle
1743 :
1744 44594 : CALL timeset(routineN, handle)
1745 :
1746 133782 : ALLOCATE (pw_grid%g(3, ng))
1747 133782 : ALLOCATE (pw_grid%gsq(ng))
1748 133782 : ALLOCATE (pw_grid%g_hat(3, ng))
1749 :
1750 44594 : nmaps = 1
1751 44594 : IF (pw_grid%grid_span == HALFSPACE) nmaps = 2
1752 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
1753 : CALL offload_activate_chosen_device()
1754 :
1755 : length = INT(int_size*MAX(ng, 1)*MAX(nmaps, 1), KIND=C_SIZE_T)
1756 : stat = offload_malloc_pinned_mem(cptr_g_hatmap, length)
1757 : CPASSERT(stat == 0)
1758 : CALL c_f_pointer(cptr_g_hatmap, pw_grid%g_hatmap, [MAX(ng, 1), MAX(nmaps, 1)])
1759 : #else
1760 178376 : ALLOCATE (pw_grid%g_hatmap(MAX(ng, 1), MAX(nmaps, 1)))
1761 : #endif
1762 :
1763 44594 : IF (pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
1764 : ALLOCATE (pw_grid%grays(pw_grid%npts(1), &
1765 157164 : pw_grid%para%nyzray(pw_grid%para%group%mepos)))
1766 : END IF
1767 :
1768 133782 : ALLOCATE (pw_grid%mapl%pos(bounds(1, 1):bounds(2, 1)))
1769 89188 : ALLOCATE (pw_grid%mapl%neg(bounds(1, 1):bounds(2, 1)))
1770 133782 : ALLOCATE (pw_grid%mapm%pos(bounds(1, 2):bounds(2, 2)))
1771 89188 : ALLOCATE (pw_grid%mapm%neg(bounds(1, 2):bounds(2, 2)))
1772 133782 : ALLOCATE (pw_grid%mapn%pos(bounds(1, 3):bounds(2, 3)))
1773 89188 : ALLOCATE (pw_grid%mapn%neg(bounds(1, 3):bounds(2, 3)))
1774 :
1775 44594 : CALL timestop(handle)
1776 :
1777 44594 : END SUBROUTINE pw_grid_allocate
1778 :
1779 : ! **************************************************************************************************
1780 : !> \brief Sort g-vectors according to length
1781 : !> \param pw_grid ...
1782 : !> \param ref_grid ...
1783 : !> \par History
1784 : !> JGH (20-12-2000) : allocate idx, ng = SIZE ( pw_grid % gsq ) the
1785 : !> sorting is local and independent from parallelisation
1786 : !> WARNING: Global ordering depends now on the number
1787 : !> of cpus.
1788 : !> JGH (28-02-2001) : check for ordering against reference grid
1789 : !> JGH (01-05-2001) : sort spherical cutoff grids also within shells
1790 : !> reference grids for non-spherical cutoffs
1791 : !> JGH (20-06-2001) : do not sort non-spherical grids
1792 : !> JGH (19-02-2003) : Order all grids, this makes subgrids also for
1793 : !> non-spherical cutoffs possible
1794 : !> JGH (21-02-2003) : Introduce gather array for general reference grids
1795 : !> \author apsi
1796 : !> Christopher Mundy
1797 : ! **************************************************************************************************
1798 44594 : SUBROUTINE pw_grid_sort(pw_grid, ref_grid)
1799 :
1800 : ! Argument
1801 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
1802 : TYPE(pw_grid_type), INTENT(IN), OPTIONAL :: ref_grid
1803 :
1804 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_sort'
1805 :
1806 : INTEGER :: handle, i, ig, ih, ip, is, it, ng, ngr
1807 : INTEGER, ALLOCATABLE, DIMENSION(:) :: idx
1808 44594 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: int_tmp
1809 : LOGICAL :: g_found
1810 : REAL(KIND=dp) :: gig, gigr
1811 44594 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: real_tmp
1812 :
1813 44594 : CALL timeset(routineN, handle)
1814 :
1815 44594 : ng = SIZE(pw_grid%gsq)
1816 133782 : ALLOCATE (idx(ng))
1817 :
1818 : ! grids are (locally) ordered by length of G-vectors
1819 44594 : CALL sort(pw_grid%gsq, ng, idx)
1820 : ! within shells order wrt x,y,z
1821 44594 : CALL sort_shells(pw_grid%gsq, pw_grid%g_hat, idx)
1822 :
1823 133782 : ALLOCATE (real_tmp(3, ng))
1824 1131568758 : DO i = 1, ng
1825 1131524164 : real_tmp(1, i) = pw_grid%g(1, idx(i))
1826 1131524164 : real_tmp(2, i) = pw_grid%g(2, idx(i))
1827 1131568758 : real_tmp(3, i) = pw_grid%g(3, idx(i))
1828 : END DO
1829 1131568758 : DO i = 1, ng
1830 1131524164 : pw_grid%g(1, i) = real_tmp(1, i)
1831 1131524164 : pw_grid%g(2, i) = real_tmp(2, i)
1832 1131568758 : pw_grid%g(3, i) = real_tmp(3, i)
1833 : END DO
1834 44594 : DEALLOCATE (real_tmp)
1835 :
1836 133782 : ALLOCATE (int_tmp(3, ng))
1837 1131568758 : DO i = 1, ng
1838 1131524164 : int_tmp(1, i) = pw_grid%g_hat(1, idx(i))
1839 1131524164 : int_tmp(2, i) = pw_grid%g_hat(2, idx(i))
1840 1131568758 : int_tmp(3, i) = pw_grid%g_hat(3, idx(i))
1841 : END DO
1842 1131568758 : DO i = 1, ng
1843 1131524164 : pw_grid%g_hat(1, i) = int_tmp(1, i)
1844 1131524164 : pw_grid%g_hat(2, i) = int_tmp(2, i)
1845 1131568758 : pw_grid%g_hat(3, i) = int_tmp(3, i)
1846 : END DO
1847 44594 : DEALLOCATE (int_tmp)
1848 :
1849 44594 : DEALLOCATE (idx)
1850 :
1851 : ! check if ordering is compatible to reference grid
1852 44594 : IF (PRESENT(ref_grid)) THEN
1853 27156 : ngr = SIZE(ref_grid%gsq)
1854 27156 : ngr = MIN(ng, ngr)
1855 27156 : IF (pw_grid%spherical) THEN
1856 0 : IF (.NOT. ALL(pw_grid%g_hat(1:3, 1:ngr) &
1857 : == ref_grid%g_hat(1:3, 1:ngr))) THEN
1858 0 : CPABORT("G space sorting not compatible")
1859 : END IF
1860 : ELSE
1861 81468 : ALLOCATE (pw_grid%gidx(1:ngr))
1862 257950330 : pw_grid%gidx = 0
1863 : ! first try as many trivial associations as possible
1864 27156 : it = 0
1865 137118780 : DO ig = 1, ngr
1866 548397222 : IF (.NOT. ALL(pw_grid%g_hat(1:3, ig) &
1867 : == ref_grid%g_hat(1:3, ig))) EXIT
1868 137091624 : pw_grid%gidx(ig) = ig
1869 137118780 : it = ig
1870 : END DO
1871 : ! now for the ones that could not be done
1872 27156 : IF (ng == ngr) THEN
1873 : ! for the case pw_grid <= ref_grid
1874 27154 : is = it
1875 120843065 : DO ig = it + 1, ngr
1876 120815911 : gig = pw_grid%gsq(ig)
1877 120815911 : gigr = MAX(1._dp, gig)
1878 120815911 : g_found = .FALSE.
1879 608284882 : DO ih = is + 1, SIZE(ref_grid%gsq)
1880 608284882 : IF (ABS(gig - ref_grid%gsq(ih))/gigr > 1.e-12_dp) CYCLE
1881 : g_found = .TRUE.
1882 487468971 : EXIT
1883 : END DO
1884 120815911 : IF (.NOT. g_found) THEN
1885 0 : WRITE (*, "(A,I10,F20.10)") "G-vector", ig, pw_grid%gsq(ig)
1886 0 : CPABORT("G vector not found")
1887 : END IF
1888 120815911 : ip = ih - 1
1889 8453566818 : DO ih = ip + 1, SIZE(ref_grid%gsq)
1890 8453566818 : IF (ABS(gig - ref_grid%gsq(ih))/gigr > 1.e-12_dp) CYCLE
1891 8453566818 : IF (pw_grid%g_hat(1, ig) /= ref_grid%g_hat(1, ih)) CYCLE
1892 457451355 : IF (pw_grid%g_hat(2, ig) /= ref_grid%g_hat(2, ih)) CYCLE
1893 149148648 : IF (pw_grid%g_hat(3, ig) /= ref_grid%g_hat(3, ih)) CYCLE
1894 120815911 : pw_grid%gidx(ig) = ih
1895 8453566818 : EXIT
1896 : END DO
1897 120815911 : IF (pw_grid%gidx(ig) == 0) THEN
1898 0 : WRITE (*, "(A,2I10)") " G-Shell ", is + 1, ip + 1
1899 : WRITE (*, "(A,I10,3I6,F20.10)") &
1900 0 : " G-vector", ig, pw_grid%g_hat(1:3, ig), pw_grid%gsq(ig)
1901 0 : DO ih = 1, SIZE(ref_grid%gsq)
1902 0 : IF (pw_grid%g_hat(1, ig) /= ref_grid%g_hat(1, ih)) CYCLE
1903 0 : IF (pw_grid%g_hat(2, ig) /= ref_grid%g_hat(2, ih)) CYCLE
1904 0 : IF (pw_grid%g_hat(3, ig) /= ref_grid%g_hat(3, ih)) CYCLE
1905 : WRITE (*, "(A,I10,3I6,F20.10)") &
1906 0 : " G-vector", ih, ref_grid%g_hat(1:3, ih), ref_grid%gsq(ih)
1907 : END DO
1908 0 : CPABORT("G vector not found")
1909 : END IF
1910 120843065 : is = ip
1911 : END DO
1912 : ELSE
1913 : ! for the case pw_grid > ref_grid
1914 2 : is = it
1915 15641 : DO ig = it + 1, ngr
1916 15639 : gig = ref_grid%gsq(ig)
1917 15639 : gigr = MAX(1._dp, gig)
1918 15639 : g_found = .FALSE.
1919 44681 : DO ih = is + 1, ng
1920 44681 : IF (ABS(pw_grid%gsq(ih) - gig)/gigr > 1.e-12_dp) CYCLE
1921 : g_found = .TRUE.
1922 29042 : EXIT
1923 : END DO
1924 15639 : IF (.NOT. g_found) THEN
1925 0 : WRITE (*, "(A,I10,F20.10)") "G-vector", ig, ref_grid%gsq(ig)
1926 0 : CPABORT("G vector not found")
1927 : END IF
1928 15639 : ip = ih - 1
1929 534662 : DO ih = ip + 1, ng
1930 534662 : IF (ABS(pw_grid%gsq(ih) - gig)/gigr > 1.e-12_dp) CYCLE
1931 534662 : IF (pw_grid%g_hat(1, ih) /= ref_grid%g_hat(1, ig)) CYCLE
1932 52928 : IF (pw_grid%g_hat(2, ih) /= ref_grid%g_hat(2, ig)) CYCLE
1933 19626 : IF (pw_grid%g_hat(3, ih) /= ref_grid%g_hat(3, ig)) CYCLE
1934 15639 : pw_grid%gidx(ig) = ih
1935 534662 : EXIT
1936 : END DO
1937 15639 : IF (pw_grid%gidx(ig) == 0) THEN
1938 0 : WRITE (*, "(A,2I10)") " G-Shell ", is + 1, ip + 1
1939 : WRITE (*, "(A,I10,3I6,F20.10)") &
1940 0 : " G-vector", ig, ref_grid%g_hat(1:3, ig), ref_grid%gsq(ig)
1941 0 : DO ih = 1, ng
1942 0 : IF (pw_grid%g_hat(1, ih) /= ref_grid%g_hat(1, ig)) CYCLE
1943 0 : IF (pw_grid%g_hat(2, ih) /= ref_grid%g_hat(2, ig)) CYCLE
1944 0 : IF (pw_grid%g_hat(3, ih) /= ref_grid%g_hat(3, ig)) CYCLE
1945 : WRITE (*, "(A,I10,3I6,F20.10)") &
1946 0 : " G-vector", ih, pw_grid%g_hat(1:3, ih), pw_grid%gsq(ih)
1947 : END DO
1948 0 : CPABORT("G vector not found")
1949 : END IF
1950 15641 : is = ip
1951 : END DO
1952 : END IF
1953 : ! test if all g-vectors are associated
1954 257950330 : IF (ANY(pw_grid%gidx == 0)) THEN
1955 0 : CPABORT("G space sorting not compatible")
1956 : END IF
1957 : END IF
1958 : END IF
1959 :
1960 : !check if G=0 is at first position
1961 44594 : IF (pw_grid%have_g0) THEN
1962 24936 : IF (pw_grid%gsq(1) /= 0.0_dp) THEN
1963 0 : CPABORT("G = 0 not in first position")
1964 : END IF
1965 : END IF
1966 :
1967 44594 : CALL timestop(handle)
1968 :
1969 44594 : END SUBROUTINE pw_grid_sort
1970 :
1971 : ! **************************************************************************************************
1972 : !> \brief ...
1973 : !> \param gsq ...
1974 : !> \param g_hat ...
1975 : !> \param idx ...
1976 : ! **************************************************************************************************
1977 44594 : SUBROUTINE sort_shells(gsq, g_hat, idx)
1978 :
1979 : ! Argument
1980 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: gsq
1981 : INTEGER, DIMENSION(:, :), INTENT(IN) :: g_hat
1982 : INTEGER, DIMENSION(:), INTENT(INOUT) :: idx
1983 :
1984 : CHARACTER(len=*), PARAMETER :: routineN = 'sort_shells'
1985 : REAL(KIND=dp), PARAMETER :: small = 5.e-16_dp
1986 :
1987 : INTEGER :: handle, ig, ng, s1, s2
1988 : REAL(KIND=dp) :: s_begin
1989 :
1990 44594 : CALL timeset(routineN, handle)
1991 :
1992 : ! Juergs temporary hack to get the grids sorted for large (4000Ry) cutoffs.
1993 : ! might need to call lapack for machine precision.
1994 :
1995 44594 : ng = SIZE(gsq)
1996 44594 : s_begin = -1.0_dp
1997 44594 : s1 = 0
1998 44594 : s2 = 0
1999 : ig = 1
2000 1131568758 : DO ig = 1, ng
2001 1131568758 : IF (ABS(gsq(ig) - s_begin) < small) THEN
2002 1065983269 : s2 = ig
2003 : ELSE
2004 65540895 : CALL redist(g_hat, idx, s1, s2)
2005 65540895 : s_begin = gsq(ig)
2006 65540895 : s1 = ig
2007 65540895 : s2 = ig
2008 : END IF
2009 : END DO
2010 44594 : CALL redist(g_hat, idx, s1, s2)
2011 :
2012 44594 : CALL timestop(handle)
2013 :
2014 44594 : END SUBROUTINE sort_shells
2015 :
2016 : ! **************************************************************************************************
2017 : !> \brief ...
2018 : !> \param g_hat ...
2019 : !> \param idx ...
2020 : !> \param s1 ...
2021 : !> \param s2 ...
2022 : ! **************************************************************************************************
2023 65585489 : SUBROUTINE redist(g_hat, idx, s1, s2)
2024 :
2025 : ! Argument
2026 : INTEGER, DIMENSION(:, :), INTENT(IN) :: g_hat
2027 : INTEGER, DIMENSION(:), INTENT(INOUT) :: idx
2028 : INTEGER, INTENT(IN) :: s1, s2
2029 :
2030 : INTEGER :: i, ii, n1, n2, n3, ns
2031 65585489 : INTEGER, ALLOCATABLE, DIMENSION(:) :: indl
2032 65585489 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: slen
2033 :
2034 65585489 : IF (s2 <= s1) RETURN
2035 61991583 : ns = s2 - s1 + 1
2036 185974749 : ALLOCATE (indl(ns))
2037 185974749 : ALLOCATE (slen(ns))
2038 :
2039 1189966435 : DO i = s1, s2
2040 1127974852 : ii = idx(i)
2041 1127974852 : n1 = g_hat(1, ii)
2042 1127974852 : n2 = g_hat(2, ii)
2043 1127974852 : n3 = g_hat(3, ii)
2044 : slen(i - s1 + 1) = 1000.0_dp*REAL(n1, dp) + &
2045 1189966435 : REAL(n2, dp) + 0.001_dp*REAL(n3, dp)
2046 : END DO
2047 61991583 : CALL sort(slen, ns, indl)
2048 1189966435 : DO i = 1, ns
2049 1127974852 : ii = indl(i) + s1 - 1
2050 1189966435 : indl(i) = idx(ii)
2051 : END DO
2052 1189966435 : idx(s1:s2) = indl(1:ns)
2053 :
2054 61991583 : DEALLOCATE (indl)
2055 61991583 : DEALLOCATE (slen)
2056 :
2057 : END SUBROUTINE redist
2058 :
2059 : ! **************************************************************************************************
2060 : !> \brief Reorder yzq and yzp arrays for parallel FFT according to FFT mapping
2061 : !> \param pw_grid ...
2062 : !> \param yz ...
2063 : !> \par History
2064 : !> none
2065 : !> \author JGH (17-Jan-2001)
2066 : ! **************************************************************************************************
2067 44594 : SUBROUTINE pw_grid_remap(pw_grid, yz)
2068 :
2069 : ! Argument
2070 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
2071 : INTEGER, DIMENSION(:, :), INTENT(OUT) :: yz
2072 :
2073 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_remap'
2074 :
2075 : INTEGER :: gpt, handle, i, ip, is, j, m, n
2076 :
2077 44594 : IF (pw_grid%para%mode == PW_MODE_LOCAL) RETURN
2078 :
2079 39291 : CALL timeset(routineN, handle)
2080 :
2081 : ASSOCIATE (ny => pw_grid%npts(2), nz => pw_grid%npts(3), posm => pw_grid%mapm%pos, posn => pw_grid%mapn%pos, &
2082 : negm => pw_grid%mapm%neg, negn => pw_grid%mapn%neg)
2083 :
2084 42030569 : yz = 0
2085 119629575 : pw_grid%para%yzp = 0
2086 42030569 : pw_grid%para%yzq = 0
2087 :
2088 1082905344 : DO gpt = 1, SIZE(pw_grid%gsq)
2089 1082866053 : m = posm(pw_grid%g_hat(2, gpt)) + 1
2090 1082866053 : n = posn(pw_grid%g_hat(3, gpt)) + 1
2091 1082905344 : yz(m, n) = yz(m, n) + 1
2092 : END DO
2093 39291 : IF (pw_grid%grid_span == HALFSPACE) THEN
2094 18717074 : DO gpt = 1, SIZE(pw_grid%gsq)
2095 18714538 : m = negm(pw_grid%g_hat(2, gpt)) + 1
2096 18714538 : n = negn(pw_grid%g_hat(3, gpt)) + 1
2097 18717074 : yz(m, n) = yz(m, n) + 1
2098 : END DO
2099 : END IF
2100 :
2101 39291 : ip = pw_grid%para%group%mepos
2102 39291 : is = 0
2103 1098737 : DO i = 1, nz
2104 42030569 : DO j = 1, ny
2105 41991278 : IF (yz(j, i) > 0) THEN
2106 19898585 : is = is + 1
2107 19898585 : pw_grid%para%yzp(1, is, ip) = j
2108 19898585 : pw_grid%para%yzp(2, is, ip) = i
2109 19898585 : pw_grid%para%yzq(j, i) = is
2110 : END IF
2111 : END DO
2112 : END DO
2113 : END ASSOCIATE
2114 :
2115 39291 : CPASSERT(is == pw_grid%para%nyzray(ip))
2116 39291 : CALL pw_grid%para%group%sum(pw_grid%para%yzp)
2117 :
2118 39291 : CALL timestop(handle)
2119 :
2120 39291 : END SUBROUTINE pw_grid_remap
2121 :
2122 : ! **************************************************************************************************
2123 : !> \brief Recalculate the g-vectors after a change of the box
2124 : !> \param cell_hmat ...
2125 : !> \param pw_grid ...
2126 : !> \par History
2127 : !> JGH (20-12-2000) : get local grid size from definition of g.
2128 : !> Assume that gsq is allocated.
2129 : !> Local routine, no information on distribution of
2130 : !> PW required.
2131 : !> JGH (8-Mar-2001) : also update information on volume and grid spaceing
2132 : !> \author apsi
2133 : !> Christopher Mundy
2134 : ! **************************************************************************************************
2135 13880 : SUBROUTINE pw_grid_change(cell_hmat, pw_grid)
2136 : ! Argument
2137 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: cell_hmat
2138 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
2139 :
2140 : INTEGER :: gpt
2141 : REAL(KIND=dp) :: cell_deth, l, m, n
2142 : REAL(KIND=dp), DIMENSION(3, 3) :: cell_h_inv
2143 13880 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: g
2144 :
2145 13880 : cell_deth = ABS(det_3x3(cell_hmat))
2146 13880 : IF (cell_deth < 1.0E-10_dp) THEN
2147 : CALL cp_abort(__LOCATION__, &
2148 : "An invalid set of cell vectors was specified. "// &
2149 0 : "The determinant det(h) is too small")
2150 : END IF
2151 13880 : cell_h_inv = inv_3x3(cell_hmat)
2152 :
2153 13880 : g => pw_grid%g
2154 :
2155 13880 : CALL cell2grid(cell_hmat, cell_h_inv, cell_deth, pw_grid)
2156 :
2157 103426088 : DO gpt = 1, SIZE(g, 2)
2158 :
2159 103412208 : l = twopi*REAL(pw_grid%g_hat(1, gpt), KIND=dp)
2160 103412208 : m = twopi*REAL(pw_grid%g_hat(2, gpt), KIND=dp)
2161 103412208 : n = twopi*REAL(pw_grid%g_hat(3, gpt), KIND=dp)
2162 :
2163 103412208 : g(1, gpt) = l*cell_h_inv(1, 1) + m*cell_h_inv(2, 1) + n*cell_h_inv(3, 1)
2164 103412208 : g(2, gpt) = l*cell_h_inv(1, 2) + m*cell_h_inv(2, 2) + n*cell_h_inv(3, 2)
2165 103412208 : g(3, gpt) = l*cell_h_inv(1, 3) + m*cell_h_inv(2, 3) + n*cell_h_inv(3, 3)
2166 :
2167 : pw_grid%gsq(gpt) = g(1, gpt)*g(1, gpt) &
2168 : + g(2, gpt)*g(2, gpt) &
2169 103426088 : + g(3, gpt)*g(3, gpt)
2170 :
2171 : END DO
2172 :
2173 13880 : END SUBROUTINE pw_grid_change
2174 :
2175 : ! **************************************************************************************************
2176 : !> \brief retains the given pw grid
2177 : !> \param pw_grid the pw grid to retain
2178 : !> \par History
2179 : !> 04.2003 created [fawzi]
2180 : !> \author fawzi
2181 : !> \note
2182 : !> see doc/ReferenceCounting.html
2183 : ! **************************************************************************************************
2184 9457361 : SUBROUTINE pw_grid_retain(pw_grid)
2185 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
2186 :
2187 9457361 : CPASSERT(pw_grid%ref_count > 0)
2188 9457361 : pw_grid%ref_count = pw_grid%ref_count + 1
2189 9457361 : END SUBROUTINE pw_grid_retain
2190 :
2191 : ! **************************************************************************************************
2192 : !> \brief releases the given pw grid
2193 : !> \param pw_grid the pw grid to release
2194 : !> \par History
2195 : !> 04.2003 created [fawzi]
2196 : !> \author fawzi
2197 : !> \note
2198 : !> see doc/ReferenceCounting.html
2199 : ! **************************************************************************************************
2200 19453922 : SUBROUTINE pw_grid_release(pw_grid)
2201 :
2202 : TYPE(pw_grid_type), POINTER :: pw_grid
2203 :
2204 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2205 : INTEGER, POINTER :: dummy_ptr
2206 : INTEGER :: stat
2207 : #endif
2208 :
2209 19453922 : IF (ASSOCIATED(pw_grid)) THEN
2210 9595709 : CPASSERT(pw_grid%ref_count > 0)
2211 9595709 : pw_grid%ref_count = pw_grid%ref_count - 1
2212 9595709 : IF (pw_grid%ref_count == 0) THEN
2213 138363 : IF (ASSOCIATED(pw_grid%gidx)) THEN
2214 27152 : DEALLOCATE (pw_grid%gidx)
2215 : END IF
2216 138363 : IF (ASSOCIATED(pw_grid%g)) THEN
2217 44589 : DEALLOCATE (pw_grid%g)
2218 : END IF
2219 138363 : IF (ASSOCIATED(pw_grid%gsq)) THEN
2220 44589 : DEALLOCATE (pw_grid%gsq)
2221 : END IF
2222 138363 : IF (ALLOCATED(pw_grid%g_hat)) THEN
2223 44589 : DEALLOCATE (pw_grid%g_hat)
2224 : END IF
2225 138363 : IF (ASSOCIATED(pw_grid%g_hatmap)) THEN
2226 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2227 : dummy_ptr => pw_grid%g_hatmap(1, 1)
2228 : stat = offload_free_pinned_mem(c_loc(dummy_ptr))
2229 : CPASSERT(stat == 0)
2230 : #else
2231 44589 : DEALLOCATE (pw_grid%g_hatmap)
2232 : #endif
2233 : END IF
2234 138363 : IF (ASSOCIATED(pw_grid%grays)) THEN
2235 39286 : DEALLOCATE (pw_grid%grays)
2236 : END IF
2237 138363 : IF (ALLOCATED(pw_grid%mapl%pos)) THEN
2238 44589 : DEALLOCATE (pw_grid%mapl%pos)
2239 : END IF
2240 138363 : IF (ALLOCATED(pw_grid%mapm%pos)) THEN
2241 44589 : DEALLOCATE (pw_grid%mapm%pos)
2242 : END IF
2243 138363 : IF (ALLOCATED(pw_grid%mapn%pos)) THEN
2244 44589 : DEALLOCATE (pw_grid%mapn%pos)
2245 : END IF
2246 138363 : IF (ALLOCATED(pw_grid%mapl%neg)) THEN
2247 44589 : DEALLOCATE (pw_grid%mapl%neg)
2248 : END IF
2249 138363 : IF (ALLOCATED(pw_grid%mapm%neg)) THEN
2250 44589 : DEALLOCATE (pw_grid%mapm%neg)
2251 : END IF
2252 138363 : IF (ALLOCATED(pw_grid%mapn%neg)) THEN
2253 44589 : DEALLOCATE (pw_grid%mapn%neg)
2254 : END IF
2255 138363 : IF (ALLOCATED(pw_grid%para%bo)) THEN
2256 44589 : DEALLOCATE (pw_grid%para%bo)
2257 : END IF
2258 138363 : IF (pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
2259 39294 : IF (ALLOCATED(pw_grid%para%yzp)) THEN
2260 39286 : DEALLOCATE (pw_grid%para%yzp)
2261 : END IF
2262 39294 : IF (ALLOCATED(pw_grid%para%yzq)) THEN
2263 39286 : DEALLOCATE (pw_grid%para%yzq)
2264 : END IF
2265 39294 : IF (ALLOCATED(pw_grid%para%nyzray)) THEN
2266 39286 : DEALLOCATE (pw_grid%para%nyzray)
2267 : END IF
2268 : END IF
2269 : ! also release groups
2270 138363 : CALL pw_grid%para%group%free()
2271 138363 : IF (ALLOCATED(pw_grid%para%pos_of_x)) THEN
2272 45601 : DEALLOCATE (pw_grid%para%pos_of_x)
2273 : END IF
2274 :
2275 138363 : IF (ASSOCIATED(pw_grid)) THEN
2276 138363 : DEALLOCATE (pw_grid)
2277 : END IF
2278 : END IF
2279 : END IF
2280 19453922 : NULLIFY (pw_grid)
2281 19453922 : END SUBROUTINE pw_grid_release
2282 :
2283 : END MODULE pw_grids
|
