Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 60268 : 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 60268 : CPASSERT(.NOT. ASSOCIATED(pw_grid))
104 3555812 : ALLOCATE (pw_grid)
105 602680 : pw_grid%bounds = bounds
106 602680 : pw_grid%bounds_local = bounds
107 241072 : pw_grid%npts = bounds(2, :) - bounds(1, :) + 1
108 241072 : pw_grid%npts_local = pw_grid%npts
109 241072 : pw_grid%ngpts = PRODUCT(INT(pw_grid%npts, KIND=int_8))
110 60268 : pw_grid%ngpts_cut = pw_grid%ngpts
111 241072 : pw_grid%ngpts_local = PRODUCT(pw_grid%npts)
112 60268 : pw_grid%ngpts_cut_local = pw_grid%ngpts_local
113 60268 : pw_grid%grid_span = FULLSPACE
114 60268 : pw_grid%para%mode = PW_MODE_LOCAL
115 60268 : pw_grid%reference = 0
116 60268 : pw_grid%ref_count = 1
117 60268 : NULLIFY (pw_grid%g)
118 60268 : NULLIFY (pw_grid%gsq)
119 60268 : NULLIFY (pw_grid%g_hatmap)
120 60268 : NULLIFY (pw_grid%gidx)
121 60268 : NULLIFY (pw_grid%grays)
122 :
123 : ! assign a unique tag to this grid
124 60268 : grid_tag = grid_tag + 1
125 60268 : pw_grid%id_nr = grid_tag
126 :
127 : ! parallel info
128 180804 : rs_dims = 1
129 60268 : CALL pw_grid%para%group%create(mp_comm_self, 2, rs_dims)
130 60268 : IF (pw_grid%para%group%num_pe > 1) THEN
131 0 : pw_grid%para%mode = PW_MODE_DISTRIBUTED
132 : ELSE
133 60268 : pw_grid%para%mode = PW_MODE_LOCAL
134 : END IF
135 :
136 60268 : 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 7396894 : 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 7396894 : 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 18 : equal = .FALSE.
160 : END IF
161 :
162 7396894 : 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 64702 : 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 64702 : CPASSERT(pw_grid%ref_count > 0)
198 :
199 64702 : IF (PRESENT(id_nr)) id_nr = pw_grid%id_nr
200 64702 : IF (PRESENT(mode)) mode = pw_grid%para%mode
201 64702 : IF (PRESENT(vol)) vol = pw_grid%vol
202 64702 : IF (PRESENT(dvol)) dvol = pw_grid%dvol
203 323334 : IF (PRESENT(npts)) npts(1:3) = pw_grid%npts(1:3)
204 64702 : IF (PRESENT(ngpts)) ngpts = pw_grid%ngpts
205 64702 : IF (PRESENT(ngpts_cut)) ngpts_cut = pw_grid%ngpts_cut
206 64838 : IF (PRESENT(dr)) dr = pw_grid%dr
207 64702 : IF (PRESENT(cutoff)) cutoff = pw_grid%cutoff
208 64702 : IF (PRESENT(orthorhombic)) orthorhombic = pw_grid%orthorhombic
209 64702 : IF (PRESENT(gvectors)) gvectors => pw_grid%g
210 64702 : IF (PRESENT(gsquare)) gsquare => pw_grid%gsq
211 :
212 64702 : 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 67894 : 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 67894 : CPASSERT(pw_grid%ref_count > 0)
236 :
237 67894 : IF (PRESENT(grid_span)) THEN
238 33349 : pw_grid%grid_span = grid_span
239 : END IF
240 67894 : 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 67894 : ELSE IF (PRESENT(bounds)) THEN
245 11190 : pw_grid%bounds = bounds
246 4476 : pw_grid%npts = bounds(2, :) - bounds(1, :) + 1
247 66775 : ELSE IF (PRESENT(npts)) THEN
248 133704 : pw_grid%npts = npts
249 334260 : pw_grid%bounds = pw_grid_bounds_from_n(npts)
250 : END IF
251 67894 : IF (PRESENT(cutoff)) THEN
252 34545 : pw_grid%cutoff = cutoff
253 34545 : IF (PRESENT(spherical)) THEN
254 34545 : pw_grid%spherical = spherical
255 : ELSE
256 0 : pw_grid%spherical = .FALSE.
257 : END IF
258 : END IF
259 :
260 67894 : 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 34545 : 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 34545 : CALL timeset(routineN, handle)
313 :
314 34545 : CPASSERT(.NOT. ASSOCIATED(pw_grid))
315 1796340 : ALLOCATE (pw_grid)
316 345450 : pw_grid%bounds = 0
317 34545 : pw_grid%cutoff = 0.0_dp
318 34545 : pw_grid%grid_span = FULLSPACE
319 : pw_grid%para%mode = PW_MODE_LOCAL
320 : pw_grid%reference = 0
321 34545 : pw_grid%ref_count = 1
322 34545 : NULLIFY (pw_grid%g)
323 34545 : NULLIFY (pw_grid%gsq)
324 34545 : NULLIFY (pw_grid%g_hatmap)
325 34545 : NULLIFY (pw_grid%gidx)
326 34545 : NULLIFY (pw_grid%grays)
327 :
328 : ! assign a unique tag to this grid
329 34545 : grid_tag = grid_tag + 1
330 34545 : pw_grid%id_nr = grid_tag
331 :
332 : ! parallel info
333 34545 : IF (mp_comm%num_pe > 1) THEN
334 29328 : pw_grid%para%mode = PW_MODE_DISTRIBUTED
335 : ELSE
336 : pw_grid%para%mode = PW_MODE_LOCAL
337 : END IF
338 :
339 34545 : cell_deth = ABS(det_3x3(cell_hmat))
340 34545 : 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 34545 : cell_h_inv = inv_3x3(cell_hmat)
346 :
347 34545 : IF (PRESENT(grid_span)) THEN
348 33349 : CALL set_pw_grid_info(pw_grid, grid_span=grid_span)
349 : END IF
350 :
351 34545 : IF (PRESENT(spherical)) THEN
352 34331 : my_spherical = spherical
353 : ELSE
354 214 : my_spherical = .FALSE.
355 : END IF
356 :
357 34545 : IF (PRESENT(odd)) THEN
358 29936 : my_odd = odd
359 : ELSE
360 4609 : my_odd = .FALSE.
361 : END IF
362 :
363 34545 : IF (PRESENT(fft_usage)) THEN
364 34419 : my_fft_usage = fft_usage
365 : ELSE
366 126 : my_fft_usage = .FALSE.
367 : END IF
368 :
369 34545 : IF (PRESENT(ncommensurate)) THEN
370 29876 : my_ncommensurate = ncommensurate
371 29876 : IF (PRESENT(icommensurate)) THEN
372 29876 : my_icommensurate = icommensurate
373 : ELSE
374 0 : my_icommensurate = MIN(1, ncommensurate)
375 : END IF
376 : ELSE
377 4669 : my_ncommensurate = 0
378 4669 : my_icommensurate = 1
379 : END IF
380 :
381 34545 : IF (PRESENT(bounds)) THEN
382 1119 : IF (PRESENT(cutoff)) THEN
383 : CALL set_pw_grid_info(pw_grid, bounds=bounds, cutoff=cutoff, &
384 172 : spherical=my_spherical)
385 : ELSE
386 3788 : n = bounds(2, :) - bounds(1, :) + 1
387 947 : my_cutoff = pw_find_cutoff(n, cell_h_inv)
388 947 : my_cutoff = 0.5_dp*my_cutoff*my_cutoff
389 : CALL set_pw_grid_info(pw_grid, bounds=bounds, cutoff=my_cutoff, &
390 947 : spherical=my_spherical)
391 : END IF
392 33426 : ELSE IF (PRESENT(npts)) THEN
393 2398 : n = npts
394 2398 : IF (PRESENT(cutoff)) THEN
395 22 : my_cutoff = cutoff
396 : ELSE
397 2376 : my_cutoff = pw_find_cutoff(npts, cell_h_inv)
398 2376 : my_cutoff = 0.5_dp*my_cutoff*my_cutoff
399 : END IF
400 2398 : 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 9592 : ref_grid=ref_grid, n_orig=n)
405 : END IF
406 : CALL set_pw_grid_info(pw_grid, npts=n, cutoff=my_cutoff, &
407 2398 : spherical=my_spherical)
408 31028 : 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 31028 : ref_grid=ref_grid)
413 : CALL set_pw_grid_info(pw_grid, npts=n, cutoff=cutoff, &
414 31028 : 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 34545 : blocked=blocked, ref_grid=ref_grid, rs_dims=rs_dims, iounit=iounit)
421 :
422 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
423 : CALL pw_grid_create_ghatmap(pw_grid)
424 : #endif
425 :
426 34545 : CALL timestop(handle)
427 :
428 34545 : END SUBROUTINE pw_grid_create_extended
429 :
430 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
431 : ! **************************************************************************************************
432 : !> \brief sets up a combined index for CUDA gather and scatter
433 : !> \param pw_grid ...
434 : !> \author Gloess Andreas (xx-Dec-2012)
435 : ! **************************************************************************************************
436 : SUBROUTINE pw_grid_create_ghatmap(pw_grid)
437 :
438 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
439 :
440 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_create_ghatmap'
441 :
442 : INTEGER :: gpt, handle, l, m, mn, n
443 :
444 : CALL timeset(routineN, handle)
445 :
446 : ! some checks
447 : CPASSERT(pw_grid%ref_count > 0)
448 :
449 : ! mapping of map_x( g_hat(i,j)) to g_hatmap
450 : ! the second index is for switching from gather(1) to scatter(2)
451 : ASSOCIATE (g_hat => pw_grid%g_hat, g_hatmap => pw_grid%g_hatmap, pmapl => pw_grid%mapl%pos, &
452 : pmapm => pw_grid%mapm%pos, pmapn => pw_grid%mapn%pos, nmapl => pw_grid%mapl%neg, &
453 : nmapm => pw_grid%mapm%neg, nmapn => pw_grid%mapn%neg, ngpts => SIZE(pw_grid%gsq), &
454 : npts => pw_grid%npts, yzq => pw_grid%para%yzq)
455 : ! initialize map array to minus one, to guarantee memory
456 : ! range checking errors in CUDA part (just to be sure)
457 : g_hatmap(:, :) = -1
458 : IF (pw_grid%para%mode /= PW_MODE_DISTRIBUTED) THEN
459 : DO gpt = 1, ngpts
460 : l = pmapl(g_hat(1, gpt))
461 : m = pmapm(g_hat(2, gpt))
462 : n = pmapn(g_hat(3, gpt))
463 : !ATTENTION: C-mapping [start-index=0] !!!!
464 : !ATTENTION: potential integer overflow !!!!
465 : g_hatmap(gpt, 1) = l + npts(1)*(m + npts(2)*n)
466 : END DO
467 : IF (pw_grid%grid_span == HALFSPACE) THEN
468 : DO gpt = 1, ngpts
469 : l = nmapl(g_hat(1, gpt))
470 : m = nmapm(g_hat(2, gpt))
471 : n = nmapn(g_hat(3, gpt))
472 : !ATTENTION: C-mapping [start-index=0] !!!!
473 : !ATTENTION: potential integer overflow !!!!
474 : g_hatmap(gpt, 2) = l + npts(1)*(m + npts(2)*n)
475 : END DO
476 : END IF
477 : ELSE
478 : DO gpt = 1, ngpts
479 : l = pmapl(g_hat(1, gpt))
480 : m = pmapm(g_hat(2, gpt)) + 1
481 : n = pmapn(g_hat(3, gpt)) + 1
482 : !ATTENTION: C-mapping [start-index=0] !!!!
483 : !ATTENTION: potential integer overflow !!!!
484 : mn = yzq(m, n) - 1
485 : g_hatmap(gpt, 1) = l + npts(1)*mn
486 : END DO
487 : IF (pw_grid%grid_span == HALFSPACE) THEN
488 : DO gpt = 1, ngpts
489 : l = nmapl(g_hat(1, gpt))
490 : m = nmapm(g_hat(2, gpt)) + 1
491 : n = nmapn(g_hat(3, gpt)) + 1
492 : !ATTENTION: C-mapping [start-index=0] !!!!
493 : !ATTENTION: potential integer overflow !!!!
494 : mn = yzq(m, n) - 1
495 : g_hatmap(gpt, 2) = l + npts(1)*mn
496 : END DO
497 : END IF
498 : END IF
499 : END ASSOCIATE
500 :
501 : CALL timestop(handle)
502 :
503 : END SUBROUTINE pw_grid_create_ghatmap
504 : #endif
505 :
506 : ! **************************************************************************************************
507 : !> \brief sets up a pw_grid, needs valid bounds as input, it is up to you to
508 : !> make sure of it using pw_grid_bounds_from_n
509 : !> \param cell_hmat ...
510 : !> \param cell_h_inv ...
511 : !> \param cell_deth ...
512 : !> \param pw_grid ...
513 : !> \param mp_comm ...
514 : !> \param bounds_local ...
515 : !> \param blocked ...
516 : !> \param ref_grid ...
517 : !> \param rs_dims ...
518 : !> \param iounit ...
519 : !> \par History
520 : !> JGH (20-Dec-2000) : Adapted for parallel use
521 : !> JGH (28-Feb-2001) : New optional argument fft_usage
522 : !> JGH (21-Mar-2001) : Reference grid code
523 : !> JGH (21-Mar-2001) : New optional argument symm_usage
524 : !> JGH (22-Mar-2001) : Simplify group assignment (mp_comm_dup)
525 : !> JGH (21-May-2002) : Remove orthorhombic keyword (code is fast enough)
526 : !> JGH (19-Feb-2003) : Negative cutoff can be used for non-spheric grids
527 : !> Joost VandeVondele (Feb-2004) : optionally generate pw grids that are commensurate in rs
528 : !> JGH (18-Dec-2007) : Refactoring
529 : !> \author fawzi
530 : !> \note
531 : !> this is the function that should be used in the future
532 : ! **************************************************************************************************
533 34545 : SUBROUTINE pw_grid_setup_internal(cell_hmat, cell_h_inv, cell_deth, pw_grid, mp_comm, bounds_local, &
534 : blocked, ref_grid, rs_dims, iounit)
535 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: cell_hmat, cell_h_inv
536 : REAL(KIND=dp), INTENT(IN) :: cell_deth
537 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
538 :
539 : CLASS(mp_comm_type), INTENT(IN) :: mp_comm
540 : INTEGER, DIMENSION(2, 3), INTENT(IN), OPTIONAL :: bounds_local
541 : INTEGER, INTENT(in), OPTIONAL :: blocked
542 : TYPE(pw_grid_type), INTENT(in), OPTIONAL :: ref_grid
543 : INTEGER, DIMENSION(2), INTENT(in), OPTIONAL :: rs_dims
544 : INTEGER, INTENT(in), OPTIONAL :: iounit
545 :
546 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_setup_internal'
547 :
548 : INTEGER :: handle, n(3)
549 34545 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: yz_mask
550 : REAL(KIND=dp) :: ecut
551 :
552 : !------------------------------------------------------------------------------
553 :
554 34545 : CALL timeset(routineN, handle)
555 :
556 34545 : CPASSERT(pw_grid%ref_count > 0)
557 :
558 : ! set pointer to possible reference grid
559 34545 : IF (PRESENT(ref_grid)) THEN
560 21076 : pw_grid%reference = ref_grid%id_nr
561 : END IF
562 :
563 34545 : IF (pw_grid%spherical) THEN
564 3237 : ecut = pw_grid%cutoff
565 : ELSE
566 31308 : ecut = 1.e10_dp
567 : END IF
568 :
569 138180 : n(:) = pw_grid%npts(:)
570 :
571 : ! Find the number of grid points
572 : ! yz_mask counts the number of g-vectors orthogonal to the yz plane
573 : ! the indices in yz_mask are from -n/2 .. n/2 shifted by n/2 + 1
574 : ! these are not mapped indices !
575 138180 : ALLOCATE (yz_mask(n(2), n(3)))
576 34545 : CALL pw_grid_count(cell_h_inv, pw_grid, mp_comm, ecut, yz_mask)
577 :
578 : ! Check if reference grid is compatible
579 34545 : IF (PRESENT(ref_grid)) THEN
580 21076 : CPASSERT(pw_grid%para%mode == ref_grid%para%mode)
581 21076 : CPASSERT(pw_grid%grid_span == ref_grid%grid_span)
582 21076 : CPASSERT(pw_grid%spherical .EQV. ref_grid%spherical)
583 : END IF
584 :
585 : ! Distribute grid
586 : CALL pw_grid_distribute(pw_grid, mp_comm, yz_mask, bounds_local=bounds_local, ref_grid=ref_grid, blocked=blocked, &
587 34545 : rs_dims=rs_dims)
588 :
589 : ! Allocate the grid fields
590 : CALL pw_grid_allocate(pw_grid, pw_grid%ngpts_cut_local, &
591 34545 : pw_grid%bounds)
592 :
593 : ! Fill in the grid structure
594 34545 : CALL pw_grid_assign(cell_h_inv, pw_grid, ecut)
595 :
596 : ! Sort g vector wrt length (only local for each processor)
597 34545 : CALL pw_grid_sort(pw_grid, ref_grid)
598 :
599 34545 : CALL pw_grid_remap(pw_grid, yz_mask)
600 :
601 34545 : DEALLOCATE (yz_mask)
602 :
603 34545 : CALL cell2grid(cell_hmat, cell_h_inv, cell_deth, pw_grid)
604 : !
605 : ! Output: All the information of this grid type
606 : !
607 :
608 34545 : IF (PRESENT(iounit)) THEN
609 34523 : CALL pw_grid_print(pw_grid, iounit)
610 : END IF
611 :
612 34545 : CALL timestop(handle)
613 :
614 34545 : END SUBROUTINE pw_grid_setup_internal
615 :
616 : ! **************************************************************************************************
617 : !> \brief Helper routine used by pw_grid_setup_internal and pw_grid_change
618 : !> \param cell_hmat ...
619 : !> \param cell_h_inv ...
620 : !> \param cell_deth ...
621 : !> \param pw_grid ...
622 : !> \par History moved common code into new subroutine
623 : !> \author Ole Schuett
624 : ! **************************************************************************************************
625 45961 : SUBROUTINE cell2grid(cell_hmat, cell_h_inv, cell_deth, pw_grid)
626 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: cell_hmat, cell_h_inv
627 : REAL(KIND=dp), INTENT(IN) :: cell_deth
628 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
629 :
630 45961 : pw_grid%vol = ABS(cell_deth)
631 45961 : pw_grid%dvol = pw_grid%vol/REAL(pw_grid%ngpts, KIND=dp)
632 : pw_grid%dr(1) = SQRT(SUM(cell_hmat(:, 1)**2)) &
633 183844 : /REAL(pw_grid%npts(1), KIND=dp)
634 : pw_grid%dr(2) = SQRT(SUM(cell_hmat(:, 2)**2)) &
635 183844 : /REAL(pw_grid%npts(2), KIND=dp)
636 : pw_grid%dr(3) = SQRT(SUM(cell_hmat(:, 3)**2)) &
637 183844 : /REAL(pw_grid%npts(3), KIND=dp)
638 183844 : pw_grid%dh(:, 1) = cell_hmat(:, 1)/REAL(pw_grid%npts(1), KIND=dp)
639 183844 : pw_grid%dh(:, 2) = cell_hmat(:, 2)/REAL(pw_grid%npts(2), KIND=dp)
640 183844 : pw_grid%dh(:, 3) = cell_hmat(:, 3)/REAL(pw_grid%npts(3), KIND=dp)
641 183844 : pw_grid%dh_inv(1, :) = cell_h_inv(1, :)*REAL(pw_grid%npts(1), KIND=dp)
642 183844 : pw_grid%dh_inv(2, :) = cell_h_inv(2, :)*REAL(pw_grid%npts(2), KIND=dp)
643 183844 : pw_grid%dh_inv(3, :) = cell_h_inv(3, :)*REAL(pw_grid%npts(3), KIND=dp)
644 :
645 : IF ((cell_hmat(1, 2) == 0.0_dp) .AND. (cell_hmat(1, 3) == 0.0_dp) .AND. &
646 : (cell_hmat(2, 1) == 0.0_dp) .AND. (cell_hmat(2, 3) == 0.0_dp) .AND. &
647 45961 : (cell_hmat(3, 1) == 0.0_dp) .AND. (cell_hmat(3, 2) == 0.0_dp)) THEN
648 37681 : pw_grid%orthorhombic = .TRUE.
649 : ELSE
650 8280 : pw_grid%orthorhombic = .FALSE.
651 : END IF
652 45961 : END SUBROUTINE cell2grid
653 :
654 : ! **************************************************************************************************
655 : !> \brief Output of information on pw_grid
656 : !> \param pw_grid ...
657 : !> \param info ...
658 : !> \author JGH[18-05-2007] from earlier versions
659 : ! **************************************************************************************************
660 34523 : SUBROUTINE pw_grid_print(pw_grid, info)
661 :
662 : TYPE(pw_grid_type), INTENT(IN) :: pw_grid
663 : INTEGER, INTENT(IN) :: info
664 :
665 : INTEGER :: i
666 : INTEGER(KIND=int_8) :: n(3)
667 : REAL(KIND=dp) :: rv(3, 3)
668 :
669 : !------------------------------------------------------------------------------
670 : !
671 : ! Output: All the information of this grid type
672 : !
673 :
674 34523 : IF (pw_grid%para%mode == PW_MODE_LOCAL) THEN
675 5217 : IF (info > 0) THEN
676 : WRITE (info, '(/,A,T71,I10)') &
677 430 : " PW_GRID| Information for grid number ", pw_grid%id_nr
678 430 : IF (pw_grid%reference > 0) THEN
679 : WRITE (info, '(A,T71,I10)') &
680 108 : " PW_GRID| Number of the reference grid ", pw_grid%reference
681 : END IF
682 430 : WRITE (info, '(" PW_GRID| Cutoff [a.u.]",T71,f10.1)') pw_grid%cutoff
683 430 : IF (pw_grid%spherical) THEN
684 200 : WRITE (info, '(A,T78,A)') " PW_GRID| spherical cutoff: ", "YES"
685 200 : WRITE (info, '(A,T71,I10)') " PW_GRID| Grid points within cutoff", &
686 400 : pw_grid%ngpts_cut
687 : ELSE
688 230 : WRITE (info, '(A,T78,A)') " PW_GRID| spherical cutoff: ", " NO"
689 : END IF
690 1720 : DO i = 1, 3
691 1290 : WRITE (info, '(A,I3,T30,2I8,T62,A,T71,I10)') " PW_GRID| Bounds ", &
692 1290 : i, pw_grid%bounds(1, I), pw_grid%bounds(2, I), &
693 3010 : "Points:", pw_grid%npts(I)
694 : END DO
695 : WRITE (info, '(A,G12.4,T50,A,T67,F14.4)') &
696 430 : " PW_GRID| Volume element (a.u.^3)", &
697 860 : pw_grid%dvol, " Volume (a.u.^3) :", pw_grid%vol
698 430 : IF (pw_grid%grid_span == HALFSPACE) THEN
699 290 : WRITE (info, '(A,T72,A)') " PW_GRID| Grid span", "HALFSPACE"
700 : ELSE
701 140 : WRITE (info, '(A,T72,A)') " PW_GRID| Grid span", "FULLSPACE"
702 : END IF
703 : END IF
704 : ELSE
705 :
706 29306 : n(1) = pw_grid%ngpts_cut_local
707 29306 : n(2) = pw_grid%ngpts_local
708 29306 : CALL pw_grid%para%group%sum(n(1:2))
709 87918 : n(3) = SUM(pw_grid%para%nyzray)
710 117224 : rv(:, 1) = REAL(n, KIND=dp)/REAL(pw_grid%para%group%num_pe, KIND=dp)
711 29306 : n(1) = pw_grid%ngpts_cut_local
712 29306 : n(2) = pw_grid%ngpts_local
713 29306 : CALL pw_grid%para%group%max(n(1:2))
714 87918 : n(3) = MAXVAL(pw_grid%para%nyzray)
715 117224 : rv(:, 2) = REAL(n, KIND=dp)
716 29306 : n(1) = pw_grid%ngpts_cut_local
717 29306 : n(2) = pw_grid%ngpts_local
718 29306 : CALL pw_grid%para%group%min(n(1:2))
719 87918 : n(3) = MINVAL(pw_grid%para%nyzray)
720 117224 : rv(:, 3) = REAL(n, KIND=dp)
721 :
722 29306 : IF (info > 0) THEN
723 : WRITE (info, '(/,A,T71,I10)') &
724 6560 : " PW_GRID| Information for grid number ", pw_grid%id_nr
725 6560 : IF (pw_grid%reference > 0) THEN
726 : WRITE (info, '(A,T71,I10)') &
727 4147 : " PW_GRID| Number of the reference grid ", pw_grid%reference
728 : END IF
729 : WRITE (info, '(A,T60,I10,A)') &
730 6560 : " PW_GRID| Grid distributed over ", pw_grid%para%group%num_pe, &
731 13120 : " processors"
732 : WRITE (info, '(A,T71,2I5)') &
733 6560 : " PW_GRID| Real space group dimensions ", pw_grid%para%group%num_pe_cart
734 6560 : IF (pw_grid%para%blocked) THEN
735 6 : WRITE (info, '(A,T78,A)') " PW_GRID| the grid is blocked: ", "YES"
736 : ELSE
737 6554 : WRITE (info, '(A,T78,A)') " PW_GRID| the grid is blocked: ", " NO"
738 : END IF
739 6560 : WRITE (info, '(" PW_GRID| Cutoff [a.u.]",T71,f10.1)') pw_grid%cutoff
740 6560 : IF (pw_grid%spherical) THEN
741 395 : WRITE (info, '(A,T78,A)') " PW_GRID| spherical cutoff: ", "YES"
742 395 : WRITE (info, '(A,T71,I10)') " PW_GRID| Grid points within cutoff", &
743 790 : pw_grid%ngpts_cut
744 : ELSE
745 6165 : WRITE (info, '(A,T78,A)') " PW_GRID| spherical cutoff: ", " NO"
746 : END IF
747 26240 : DO i = 1, 3
748 19680 : WRITE (info, '(A,I3,T30,2I8,T62,A,T71,I10)') " PW_GRID| Bounds ", &
749 19680 : i, pw_grid%bounds(1, I), pw_grid%bounds(2, I), &
750 45920 : "Points:", pw_grid%npts(I)
751 : END DO
752 : WRITE (info, '(A,G12.4,T50,A,T67,F14.4)') &
753 6560 : " PW_GRID| Volume element (a.u.^3)", &
754 13120 : pw_grid%dvol, " Volume (a.u.^3) :", pw_grid%vol
755 6560 : IF (pw_grid%grid_span == HALFSPACE) THEN
756 403 : WRITE (info, '(A,T72,A)') " PW_GRID| Grid span", "HALFSPACE"
757 : ELSE
758 6157 : WRITE (info, '(A,T72,A)') " PW_GRID| Grid span", "FULLSPACE"
759 : END IF
760 6560 : WRITE (info, '(A,T48,A)') " PW_GRID| Distribution", &
761 13120 : " Average Max Min"
762 6560 : WRITE (info, '(A,T45,F12.1,2I12)') " PW_GRID| G-Vectors", &
763 13120 : rv(1, 1), NINT(rv(1, 2)), NINT(rv(1, 3))
764 6560 : WRITE (info, '(A,T45,F12.1,2I12)') " PW_GRID| G-Rays ", &
765 13120 : rv(3, 1), NINT(rv(3, 2)), NINT(rv(3, 3))
766 6560 : WRITE (info, '(A,T45,F12.1,2I12)') " PW_GRID| Real Space Points", &
767 13120 : rv(2, 1), NINT(rv(2, 2)), NINT(rv(2, 3))
768 : END IF ! group head
769 : END IF ! local
770 :
771 34523 : END SUBROUTINE pw_grid_print
772 :
773 : ! **************************************************************************************************
774 : !> \brief Distribute grids in real and Fourier Space to the processors in group
775 : !> \param pw_grid ...
776 : !> \param mp_comm ...
777 : !> \param yz_mask ...
778 : !> \param bounds_local ...
779 : !> \param ref_grid ...
780 : !> \param blocked ...
781 : !> \param rs_dims ...
782 : !> \par History
783 : !> JGH (01-Mar-2001) optional reference grid
784 : !> JGH (22-May-2002) bug fix for pre_tag and HALFSPACE grids
785 : !> JGH (09-Sep-2003) reduce scaling for distribution
786 : !> \author JGH (22-12-2000)
787 : ! **************************************************************************************************
788 34545 : SUBROUTINE pw_grid_distribute(pw_grid, mp_comm, yz_mask, bounds_local, ref_grid, blocked, rs_dims)
789 :
790 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
791 :
792 : CLASS(mp_comm_type), INTENT(IN) :: mp_comm
793 : INTEGER, DIMENSION(:, :), INTENT(INOUT) :: yz_mask
794 : INTEGER, DIMENSION(2, 3), INTENT(IN), OPTIONAL :: bounds_local
795 : TYPE(pw_grid_type), INTENT(IN), OPTIONAL :: ref_grid
796 : INTEGER, INTENT(IN), OPTIONAL :: blocked
797 : INTEGER, DIMENSION(2), INTENT(in), OPTIONAL :: rs_dims
798 :
799 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_distribute'
800 :
801 : INTEGER :: blocked_local, coor(2), gmax, handle, i, &
802 : i1, i2, ip, ipl, ipp, itmp, j, l, lby, &
803 : lbz, lo(2), m, n, np, ns, nx, ny, nz, &
804 : rsd(2)
805 34545 : INTEGER, ALLOCATABLE, DIMENSION(:) :: pemap
806 34545 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: yz_index
807 34545 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: axis_dist_all
808 : INTEGER, DIMENSION(2) :: my_rs_dims
809 : INTEGER, DIMENSION(2, 3) :: axis_dist
810 : LOGICAL :: blocking
811 :
812 : !------------------------------------------------------------------------------
813 :
814 34545 : CALL timeset(routineN, handle)
815 :
816 34545 : lby = pw_grid%bounds(1, 2)
817 34545 : lbz = pw_grid%bounds(1, 3)
818 :
819 138180 : pw_grid%ngpts = PRODUCT(INT(pw_grid%npts, KIND=int_8))
820 :
821 34545 : my_rs_dims = 0
822 34545 : IF (PRESENT(rs_dims)) THEN
823 33532 : my_rs_dims = rs_dims
824 : END IF
825 :
826 34545 : IF (PRESENT(blocked)) THEN
827 31132 : blocked_local = blocked
828 : ELSE
829 : blocked_local = do_pw_grid_blocked_free
830 : END IF
831 :
832 34545 : pw_grid%para%blocked = .FALSE.
833 :
834 34545 : IF (pw_grid%para%mode == PW_MODE_LOCAL) THEN
835 :
836 5217 : pw_grid%para%ray_distribution = .FALSE.
837 :
838 52170 : pw_grid%bounds_local = pw_grid%bounds
839 20868 : pw_grid%npts_local = pw_grid%npts
840 5217 : CPASSERT(pw_grid%ngpts_cut < HUGE(pw_grid%ngpts_cut_local))
841 5217 : pw_grid%ngpts_cut_local = INT(pw_grid%ngpts_cut)
842 5217 : CPASSERT(pw_grid%ngpts < HUGE(pw_grid%ngpts_local))
843 5217 : pw_grid%ngpts_local = INT(pw_grid%ngpts)
844 15651 : my_rs_dims = 1
845 5217 : CALL pw_grid%para%group%create(mp_comm, 2, my_rs_dims)
846 :
847 5217 : ALLOCATE (pw_grid%para%bo(2, 3, 0:0, 3))
848 20868 : DO i = 1, 3
849 62604 : pw_grid%para%bo(1, 1:3, 0, i) = 1
850 67821 : pw_grid%para%bo(2, 1:3, 0, i) = pw_grid%npts(1:3)
851 : END DO
852 :
853 : ELSE
854 :
855 : !..find the real space distribution
856 29328 : nx = pw_grid%npts(1)
857 29328 : ny = pw_grid%npts(2)
858 29328 : nz = pw_grid%npts(3)
859 29328 : np = mp_comm%num_pe
860 :
861 : ! The user can specify 2 strictly positive indices => specific layout
862 : ! 1 strictly positive index => the other is fixed by the number of CPUs
863 : ! 0 strictly positive indices => fully free distribution
864 : ! if fully free or the user request can not be fulfilled, we optimize heuristically ourselves by:
865 : ! 1) nx>np -> taking a plane distribution (/np,1/)
866 : ! 2) nx<np -> taking the most square distribution
867 : ! if blocking is free:
868 : ! 1) blocked=.FALSE. for plane distributions
869 : ! 2) blocked=.TRUE. for non-plane distributions
870 41316 : IF (ANY(my_rs_dims <= 0)) THEN
871 69978 : IF (ALL(my_rs_dims <= 0)) THEN
872 23310 : my_rs_dims = [0, 0]
873 : ELSE
874 24 : IF (my_rs_dims(1) > 0) THEN
875 0 : my_rs_dims(2) = np/my_rs_dims(1)
876 : ELSE
877 24 : my_rs_dims(1) = np/my_rs_dims(2)
878 : END IF
879 : END IF
880 : END IF
881 : ! reset if the distribution can not be fulfilled
882 87984 : IF (PRODUCT(my_rs_dims) .NE. np) my_rs_dims = [0, 0]
883 : ! reset if the distribution can not be dealt with [1,np]
884 29376 : IF (ALL(my_rs_dims == [1, np])) my_rs_dims = [0, 0]
885 :
886 : ! if [0,0] now, we can optimize it ourselves
887 76004 : IF (ALL(my_rs_dims == [0, 0])) THEN
888 : ! only small grids have a chance to be 2d distributed
889 23338 : IF (nx < np) THEN
890 : ! gives the most square looking distribution
891 0 : CALL mp_dims_create(np, my_rs_dims)
892 : ! we tend to like the first index being smaller than the second
893 0 : IF (my_rs_dims(1) > my_rs_dims(2)) THEN
894 0 : itmp = my_rs_dims(1)
895 0 : my_rs_dims(1) = my_rs_dims(2)
896 0 : my_rs_dims(2) = itmp
897 : END IF
898 : ! but should avoid having the first index 1 in all cases
899 0 : IF (my_rs_dims(1) == 1) THEN
900 0 : itmp = my_rs_dims(1)
901 0 : my_rs_dims(1) = my_rs_dims(2)
902 0 : my_rs_dims(2) = itmp
903 : END IF
904 : ELSE
905 70014 : my_rs_dims = [np, 1]
906 : END IF
907 : END IF
908 :
909 : ! now fix the blocking if we have a choice
910 : SELECT CASE (blocked_local)
911 : CASE (do_pw_grid_blocked_false)
912 28 : blocking = .FALSE.
913 : CASE (do_pw_grid_blocked_true)
914 28 : blocking = .TRUE.
915 : CASE (do_pw_grid_blocked_free)
916 29934 : IF (ALL(my_rs_dims == [np, 1])) THEN
917 : blocking = .FALSE.
918 : ELSE
919 0 : blocking = .TRUE.
920 : END IF
921 : CASE DEFAULT
922 29328 : CPABORT("")
923 : END SELECT
924 :
925 : !..create group for real space distribution
926 29328 : CALL pw_grid%para%group%create(mp_comm, 2, my_rs_dims)
927 :
928 29328 : IF (PRESENT(bounds_local)) THEN
929 220 : pw_grid%bounds_local = bounds_local
930 : ELSE
931 : lo = get_limit(nx, pw_grid%para%group%num_pe_cart(1), &
932 29306 : pw_grid%para%group%mepos_cart(1))
933 87918 : pw_grid%bounds_local(:, 1) = lo + pw_grid%bounds(1, 1) - 1
934 : lo = get_limit(ny, pw_grid%para%group%num_pe_cart(2), &
935 29306 : pw_grid%para%group%mepos_cart(2))
936 87918 : pw_grid%bounds_local(:, 2) = lo + pw_grid%bounds(1, 2) - 1
937 87918 : pw_grid%bounds_local(:, 3) = pw_grid%bounds(:, 3)
938 : END IF
939 :
940 : pw_grid%npts_local(:) = pw_grid%bounds_local(2, :) &
941 117312 : - pw_grid%bounds_local(1, :) + 1
942 :
943 : !..the third distribution is needed for the second step in the FFT
944 117312 : ALLOCATE (pw_grid%para%bo(2, 3, 0:np - 1, 3))
945 87984 : rsd = pw_grid%para%group%num_pe_cart
946 :
947 29328 : IF (PRESENT(bounds_local)) THEN
948 : ! axis_dist tells what portion of 1 .. nx , 1 .. ny , 1 .. nz are in the current process
949 88 : DO i = 1, 3
950 220 : axis_dist(:, i) = bounds_local(:, i) - pw_grid%bounds(1, i) + 1
951 : END DO
952 66 : ALLOCATE (axis_dist_all(2, 3, np))
953 22 : CALL pw_grid%para%group%allgather(axis_dist, axis_dist_all)
954 66 : DO ip = 0, np - 1
955 44 : CALL pw_grid%para%group%coords(ip, coor)
956 : ! distribution xyZ
957 132 : pw_grid%para%bo(1:2, 1, ip, 1) = axis_dist_all(1:2, 1, ip + 1)
958 132 : pw_grid%para%bo(1:2, 2, ip, 1) = axis_dist_all(1:2, 2, ip + 1)
959 44 : pw_grid%para%bo(1, 3, ip, 1) = 1
960 44 : pw_grid%para%bo(2, 3, ip, 1) = nz
961 : ! distribution xYz
962 132 : pw_grid%para%bo(1:2, 1, ip, 2) = axis_dist_all(1:2, 1, ip + 1)
963 44 : pw_grid%para%bo(1, 2, ip, 2) = 1
964 44 : pw_grid%para%bo(2, 2, ip, 2) = ny
965 132 : pw_grid%para%bo(1:2, 3, ip, 2) = get_limit(nz, rsd(2), coor(2))
966 : ! distribution Xyz
967 44 : pw_grid%para%bo(1, 1, ip, 3) = 1
968 44 : pw_grid%para%bo(2, 1, ip, 3) = nx
969 132 : pw_grid%para%bo(1:2, 2, ip, 3) = get_limit(ny, rsd(1), coor(1))
970 154 : pw_grid%para%bo(1:2, 3, ip, 3) = get_limit(nz, rsd(2), coor(2))
971 : END DO
972 22 : DEALLOCATE (axis_dist_all)
973 : ELSE
974 87918 : DO ip = 0, np - 1
975 58612 : CALL pw_grid%para%group%coords(ip, coor)
976 : ! distribution xyZ
977 175836 : pw_grid%para%bo(1:2, 1, ip, 1) = get_limit(nx, rsd(1), coor(1))
978 175836 : pw_grid%para%bo(1:2, 2, ip, 1) = get_limit(ny, rsd(2), coor(2))
979 58612 : pw_grid%para%bo(1, 3, ip, 1) = 1
980 58612 : pw_grid%para%bo(2, 3, ip, 1) = nz
981 : ! distribution xYz
982 175836 : pw_grid%para%bo(1:2, 1, ip, 2) = get_limit(nx, rsd(1), coor(1))
983 58612 : pw_grid%para%bo(1, 2, ip, 2) = 1
984 58612 : pw_grid%para%bo(2, 2, ip, 2) = ny
985 175836 : pw_grid%para%bo(1:2, 3, ip, 2) = get_limit(nz, rsd(2), coor(2))
986 : ! distribution Xyz
987 58612 : pw_grid%para%bo(1, 1, ip, 3) = 1
988 58612 : pw_grid%para%bo(2, 1, ip, 3) = nx
989 175836 : pw_grid%para%bo(1:2, 2, ip, 3) = get_limit(ny, rsd(1), coor(1))
990 205142 : pw_grid%para%bo(1:2, 3, ip, 3) = get_limit(nz, rsd(2), coor(2))
991 : END DO
992 : END IF
993 : !..find the g space distribution
994 29328 : pw_grid%ngpts_cut_local = 0
995 :
996 87984 : ALLOCATE (pw_grid%para%nyzray(0:np - 1))
997 :
998 117312 : ALLOCATE (pw_grid%para%yzq(ny, nz))
999 :
1000 : IF (pw_grid%spherical .OR. pw_grid%grid_span == HALFSPACE &
1001 29328 : .OR. .NOT. blocking) THEN
1002 :
1003 29316 : pw_grid%para%ray_distribution = .TRUE.
1004 :
1005 31466428 : pw_grid%para%yzq = -1
1006 29316 : IF (PRESENT(ref_grid)) THEN
1007 : ! tag all vectors from the reference grid
1008 18028 : CALL pre_tag(pw_grid, yz_mask, ref_grid)
1009 : END IF
1010 :
1011 : ! Round Robin distribution
1012 : ! Processors 0 .. NP-1, NP-1 .. 0 get the largest remaining batch
1013 : ! of g vectors in turn
1014 :
1015 29316 : i1 = SIZE(yz_mask, 1)
1016 29316 : i2 = SIZE(yz_mask, 2)
1017 87948 : ALLOCATE (yz_index(2, i1*i2))
1018 29316 : CALL order_mask(yz_mask, yz_index)
1019 30697360 : DO i = 1, i1*i2
1020 30668044 : lo(1) = yz_index(1, i)
1021 30668044 : lo(2) = yz_index(2, i)
1022 30668044 : IF (lo(1)*lo(2) == 0) CYCLE
1023 18891560 : gmax = yz_mask(lo(1), lo(2))
1024 18891560 : IF (gmax == 0) CYCLE
1025 18891560 : yz_mask(lo(1), lo(2)) = 0
1026 18891560 : ip = MOD(i - 1, 2*np)
1027 18891560 : IF (ip > np - 1) ip = 2*np - ip - 1
1028 18891560 : IF (ip == pw_grid%para%group%mepos) THEN
1029 9445780 : pw_grid%ngpts_cut_local = pw_grid%ngpts_cut_local + gmax
1030 : END IF
1031 18891560 : pw_grid%para%yzq(lo(1), lo(2)) = ip
1032 18920876 : IF (pw_grid%grid_span == HALFSPACE) THEN
1033 1242744 : m = -lo(1) - 2*lby + 2
1034 1242744 : n = -lo(2) - 2*lbz + 2
1035 1242744 : pw_grid%para%yzq(m, n) = ip
1036 1242744 : yz_mask(m, n) = 0
1037 : END IF
1038 : END DO
1039 :
1040 29316 : DEALLOCATE (yz_index)
1041 :
1042 : ! Count the total number of rays on each processor
1043 87948 : pw_grid%para%nyzray = 0
1044 798384 : DO i = 1, nz
1045 31466428 : DO j = 1, ny
1046 30668044 : ip = pw_grid%para%yzq(j, i)
1047 30668044 : IF (ip >= 0) pw_grid%para%nyzray(ip) = &
1048 30223572 : pw_grid%para%nyzray(ip) + 1
1049 : END DO
1050 : END DO
1051 :
1052 : ! Allocate mapping array (y:z, nray, nproc)
1053 87948 : ns = MAXVAL(pw_grid%para%nyzray(0:np - 1))
1054 117264 : ALLOCATE (pw_grid%para%yzp(2, ns, 0:np - 1))
1055 :
1056 : ! Fill mapping array, recalculate nyzray for convenience
1057 87948 : pw_grid%para%nyzray = 0
1058 798384 : DO i = 1, nz
1059 31466428 : DO j = 1, ny
1060 30668044 : ip = pw_grid%para%yzq(j, i)
1061 31437112 : IF (ip >= 0) THEN
1062 : pw_grid%para%nyzray(ip) = &
1063 29454504 : pw_grid%para%nyzray(ip) + 1
1064 29454504 : ns = pw_grid%para%nyzray(ip)
1065 29454504 : pw_grid%para%yzp(1, ns, ip) = j
1066 29454504 : pw_grid%para%yzp(2, ns, ip) = i
1067 29454504 : IF (ip == pw_grid%para%group%mepos) THEN
1068 14727252 : pw_grid%para%yzq(j, i) = ns
1069 : ELSE
1070 14727252 : pw_grid%para%yzq(j, i) = -1
1071 : END IF
1072 : ELSE
1073 1213540 : pw_grid%para%yzq(j, i) = -2
1074 : END IF
1075 : END DO
1076 : END DO
1077 :
1078 117264 : pw_grid%ngpts_local = PRODUCT(pw_grid%npts_local)
1079 :
1080 : ELSE
1081 : !
1082 : ! block distribution of g vectors, we do not have a spherical cutoff
1083 : !
1084 :
1085 12 : pw_grid%para%blocked = .TRUE.
1086 12 : pw_grid%para%ray_distribution = .FALSE.
1087 :
1088 36 : DO ip = 0, np - 1
1089 : m = pw_grid%para%bo(2, 2, ip, 3) - &
1090 24 : pw_grid%para%bo(1, 2, ip, 3) + 1
1091 : n = pw_grid%para%bo(2, 3, ip, 3) - &
1092 24 : pw_grid%para%bo(1, 3, ip, 3) + 1
1093 36 : pw_grid%para%nyzray(ip) = n*m
1094 : END DO
1095 :
1096 12 : ipl = pw_grid%para%group%mepos
1097 : l = pw_grid%para%bo(2, 1, ipl, 3) - &
1098 12 : pw_grid%para%bo(1, 1, ipl, 3) + 1
1099 : m = pw_grid%para%bo(2, 2, ipl, 3) - &
1100 12 : pw_grid%para%bo(1, 2, ipl, 3) + 1
1101 : n = pw_grid%para%bo(2, 3, ipl, 3) - &
1102 12 : pw_grid%para%bo(1, 3, ipl, 3) + 1
1103 12 : pw_grid%ngpts_cut_local = l*m*n
1104 12 : pw_grid%ngpts_local = pw_grid%ngpts_cut_local
1105 :
1106 4816 : pw_grid%para%yzq = 0
1107 : ny = pw_grid%para%bo(2, 2, ipl, 3) - &
1108 12 : pw_grid%para%bo(1, 2, ipl, 3) + 1
1109 254 : DO n = pw_grid%para%bo(1, 3, ipl, 3), &
1110 12 : pw_grid%para%bo(2, 3, ipl, 3)
1111 242 : i = n - pw_grid%para%bo(1, 3, ipl, 3)
1112 2523 : DO m = pw_grid%para%bo(1, 2, ipl, 3), &
1113 254 : pw_grid%para%bo(2, 2, ipl, 3)
1114 2281 : j = m - pw_grid%para%bo(1, 2, ipl, 3) + 1
1115 2523 : pw_grid%para%yzq(m, n) = j + i*ny
1116 : END DO
1117 : END DO
1118 :
1119 : ! Allocate mapping array (y:z, nray, nproc)
1120 36 : ns = MAXVAL(pw_grid%para%nyzray(0:np - 1))
1121 48 : ALLOCATE (pw_grid%para%yzp(2, ns, 0:np - 1))
1122 14112 : pw_grid%para%yzp = 0
1123 :
1124 36 : ALLOCATE (pemap(0:np - 1))
1125 36 : pemap = 0
1126 12 : pemap(pw_grid%para%group%mepos) = pw_grid%para%group%mepos
1127 12 : CALL pw_grid%para%group%sum(pemap)
1128 :
1129 36 : DO ip = 0, np - 1
1130 24 : ipp = pemap(ip)
1131 24 : ns = 0
1132 508 : DO n = pw_grid%para%bo(1, 3, ipp, 3), &
1133 24 : pw_grid%para%bo(2, 3, ipp, 3)
1134 484 : i = n - pw_grid%bounds(1, 3) + 1
1135 5046 : DO m = pw_grid%para%bo(1, 2, ipp, 3), &
1136 508 : pw_grid%para%bo(2, 2, ipp, 3)
1137 4562 : j = m - pw_grid%bounds(1, 2) + 1
1138 4562 : ns = ns + 1
1139 4562 : pw_grid%para%yzp(1, ns, ip) = j
1140 5046 : pw_grid%para%yzp(2, ns, ip) = i
1141 : END DO
1142 : END DO
1143 36 : CPASSERT(ns == pw_grid%para%nyzray(ip))
1144 : END DO
1145 :
1146 12 : DEALLOCATE (pemap)
1147 :
1148 : END IF
1149 :
1150 : END IF
1151 :
1152 : ! pos_of_x(i) tells on which cpu pw%array(i,:,:) is located
1153 : ! should be computable in principle, without the need for communication
1154 34545 : IF (pw_grid%para%mode .EQ. PW_MODE_DISTRIBUTED) THEN
1155 87984 : ALLOCATE (pw_grid%para%pos_of_x(pw_grid%bounds(1, 1):pw_grid%bounds(2, 1)))
1156 788366 : pw_grid%para%pos_of_x = 0
1157 408847 : pw_grid%para%pos_of_x(pw_grid%bounds_local(1, 1):pw_grid%bounds_local(2, 1)) = pw_grid%para%group%mepos
1158 29328 : CALL pw_grid%para%group%sum(pw_grid%para%pos_of_x)
1159 : ELSE
1160 : ! this should not be needed
1161 15651 : ALLOCATE (pw_grid%para%pos_of_x(pw_grid%bounds(1, 1):pw_grid%bounds(2, 1)))
1162 88060 : pw_grid%para%pos_of_x = 0
1163 : END IF
1164 :
1165 34545 : CALL timestop(handle)
1166 :
1167 34545 : END SUBROUTINE pw_grid_distribute
1168 :
1169 : ! **************************************************************************************************
1170 : !> \brief ...
1171 : !> \param pw_grid ...
1172 : !> \param yz_mask ...
1173 : !> \param ref_grid ...
1174 : !> \par History
1175 : !> - Fix mapping bug for pw_grid eqv to ref_grid (21.11.2019, MK)
1176 : ! **************************************************************************************************
1177 18028 : SUBROUTINE pre_tag(pw_grid, yz_mask, ref_grid)
1178 :
1179 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
1180 : INTEGER, DIMENSION(:, :), INTENT(INOUT) :: yz_mask
1181 : TYPE(pw_grid_type), INTENT(IN) :: ref_grid
1182 :
1183 : INTEGER :: gmax, ig, ip, lby, lbz, my, mz, ny, nz, &
1184 : uby, ubz, y, yp, z, zp
1185 :
1186 18028 : ny = ref_grid%npts(2)
1187 18028 : nz = ref_grid%npts(3)
1188 18028 : lby = pw_grid%bounds(1, 2)
1189 18028 : lbz = pw_grid%bounds(1, 3)
1190 18028 : uby = pw_grid%bounds(2, 2)
1191 18028 : ubz = pw_grid%bounds(2, 3)
1192 18028 : my = SIZE(yz_mask, 1)
1193 18028 : mz = SIZE(yz_mask, 2)
1194 :
1195 : ! loop over all processors and all g vectors yz lines on this processor
1196 54084 : DO ip = 0, ref_grid%para%group%num_pe - 1
1197 47198548 : DO ig = 1, ref_grid%para%nyzray(ip)
1198 : ! go from mapped coordinates to original coordinates
1199 : ! 1, 2, ..., n-1, n -> 0, 1, ..., (n/2)-1, -(n/2), -(n/2)+1, ..., -2, -1
1200 47144464 : y = ref_grid%para%yzp(1, ig, ip) - 1
1201 47144464 : IF (y >= ny/2) y = y - ny
1202 47144464 : z = ref_grid%para%yzp(2, ig, ip) - 1
1203 47144464 : IF (z >= nz/2) z = z - nz
1204 : ! check if this is inside the realm of the new grid
1205 47144464 : IF (y < lby .OR. y > uby .OR. z < lbz .OR. z > ubz) CYCLE
1206 : ! go to shifted coordinates
1207 9322490 : y = y - lby + 1
1208 9322490 : z = z - lbz + 1
1209 : ! this tag is outside the cutoff range of the new grid
1210 9322490 : IF (pw_grid%grid_span == HALFSPACE) THEN
1211 0 : yp = -y - 2*lby + 2
1212 0 : zp = -z - 2*lbz + 2
1213 : ! if the reference grid is larger than the mirror point may be
1214 : ! outside the new grid even if the original point is inside
1215 0 : IF (yp < 1 .OR. yp > my .OR. zp < 1 .OR. zp > mz) CYCLE
1216 0 : gmax = MAX(yz_mask(y, z), yz_mask(yp, zp))
1217 0 : IF (gmax == 0) CYCLE
1218 0 : yz_mask(y, z) = 0
1219 0 : yz_mask(yp, zp) = 0
1220 0 : pw_grid%para%yzq(y, z) = ip
1221 0 : pw_grid%para%yzq(yp, zp) = ip
1222 : ELSE
1223 9322490 : gmax = yz_mask(y, z)
1224 9322490 : IF (gmax == 0) CYCLE
1225 9322490 : yz_mask(y, z) = 0
1226 9322490 : pw_grid%para%yzq(y, z) = ip
1227 : END IF
1228 9358546 : IF (ip == pw_grid%para%group%mepos) THEN
1229 4661245 : pw_grid%ngpts_cut_local = pw_grid%ngpts_cut_local + gmax
1230 : END IF
1231 : END DO
1232 : END DO
1233 :
1234 18028 : END SUBROUTINE pre_tag
1235 :
1236 : ! **************************************************************************************************
1237 : !> \brief ...
1238 : !> \param yz_mask ...
1239 : !> \param yz_index ...
1240 : ! **************************************************************************************************
1241 29316 : PURE SUBROUTINE order_mask(yz_mask, yz_index)
1242 :
1243 : INTEGER, DIMENSION(:, :), INTENT(IN) :: yz_mask
1244 : INTEGER, DIMENSION(:, :), INTENT(OUT) :: yz_index
1245 :
1246 : INTEGER :: i1, i2, ic, icount, ii, im, jc, jj
1247 :
1248 : !NB load balance
1249 : !------------------------------------------------------------------------------
1250 : !NB spiral out from origin, so that even if overall grid is full and
1251 : !NB block distributed, spherical cutoff still leads to good load
1252 : !NB balance in cp_ddapc_apply_CD
1253 :
1254 29316 : i1 = SIZE(yz_mask, 1)
1255 29316 : i2 = SIZE(yz_mask, 2)
1256 92033448 : yz_index = 0
1257 :
1258 29316 : icount = 1
1259 29316 : ic = i1/2
1260 29316 : jc = i2/2
1261 29316 : ii = ic
1262 29316 : jj = jc
1263 29316 : IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1264 29316 : IF (yz_mask(ii, jj) /= 0) THEN
1265 11288 : yz_index(1, icount) = ii
1266 11288 : yz_index(2, icount) = jj
1267 11288 : icount = icount + 1
1268 : END IF
1269 : END IF
1270 440582 : DO im = 1, MAX(ic + 1, jc + 1)
1271 411266 : ii = ic - im
1272 10541004 : DO jj = jc - im, jc + im
1273 10541004 : IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1274 7356522 : IF (yz_mask(ii, jj) /= 0) THEN
1275 4590702 : yz_index(1, icount) = ii
1276 4590702 : yz_index(2, icount) = jj
1277 4590702 : icount = icount + 1
1278 : END IF
1279 : END IF
1280 : END DO
1281 411266 : ii = ic + im
1282 10541004 : DO jj = jc - im, jc + im
1283 10541004 : IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1284 8313906 : IF (yz_mask(ii, jj) /= 0) THEN
1285 5009142 : yz_index(1, icount) = ii
1286 5009142 : yz_index(2, icount) = jj
1287 5009142 : icount = icount + 1
1288 : END IF
1289 : END IF
1290 : END DO
1291 411266 : jj = jc - im
1292 9718472 : DO ii = ic - im + 1, ic + im - 1
1293 9718472 : IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1294 7026084 : IF (yz_mask(ii, jj) /= 0) THEN
1295 4721078 : yz_index(1, icount) = ii
1296 4721078 : yz_index(2, icount) = jj
1297 4721078 : icount = icount + 1
1298 : END IF
1299 : END IF
1300 : END DO
1301 9718472 : jj = jc + im
1302 9747788 : DO ii = ic - im + 1, ic + im - 1
1303 9718472 : IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1304 7942216 : IF (yz_mask(ii, jj) /= 0) THEN
1305 4559350 : yz_index(1, icount) = ii
1306 4559350 : yz_index(2, icount) = jj
1307 4559350 : icount = icount + 1
1308 : END IF
1309 : END IF
1310 : END DO
1311 : END DO
1312 :
1313 29316 : END SUBROUTINE order_mask
1314 : ! **************************************************************************************************
1315 : !> \brief compute the length of g vectors
1316 : !> \param h_inv ...
1317 : !> \param length_x ...
1318 : !> \param length_y ...
1319 : !> \param length_z ...
1320 : !> \param length ...
1321 : !> \param l ...
1322 : !> \param m ...
1323 : !> \param n ...
1324 : ! **************************************************************************************************
1325 1695321345 : PURE SUBROUTINE pw_vec_length(h_inv, length_x, length_y, length_z, length, l, m, n)
1326 :
1327 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: h_inv
1328 : REAL(KIND=dp), INTENT(OUT) :: length_x, length_y, length_z, length
1329 : INTEGER, INTENT(IN) :: l, m, n
1330 :
1331 : length_x &
1332 : = REAL(l, dp)*h_inv(1, 1) &
1333 : + REAL(m, dp)*h_inv(2, 1) &
1334 1695321345 : + REAL(n, dp)*h_inv(3, 1)
1335 : length_y &
1336 : = REAL(l, dp)*h_inv(1, 2) &
1337 : + REAL(m, dp)*h_inv(2, 2) &
1338 1695321345 : + REAL(n, dp)*h_inv(3, 2)
1339 : length_z &
1340 : = REAL(l, dp)*h_inv(1, 3) &
1341 : + REAL(m, dp)*h_inv(2, 3) &
1342 1695321345 : + REAL(n, dp)*h_inv(3, 3)
1343 :
1344 : ! enforce strict zero-ness in this case (compiler optimization)
1345 1695321345 : IF (l == 0 .AND. m == 0 .AND. n == 0) THEN
1346 39762 : length_x = 0
1347 39762 : length_y = 0
1348 39762 : length_z = 0
1349 : END IF
1350 :
1351 1695321345 : length_x = length_x*twopi
1352 1695321345 : length_y = length_y*twopi
1353 1695321345 : length_z = length_z*twopi
1354 :
1355 1695321345 : length = length_x**2 + length_y**2 + length_z**2
1356 :
1357 1695321345 : END SUBROUTINE
1358 :
1359 : ! **************************************************************************************************
1360 : !> \brief Count total number of g vectors
1361 : !> \param h_inv ...
1362 : !> \param pw_grid ...
1363 : !> \param mp_comm ...
1364 : !> \param cutoff ...
1365 : !> \param yz_mask ...
1366 : !> \par History
1367 : !> JGH (22-12-2000) : Adapted for parallel use
1368 : !> \author apsi
1369 : !> Christopher Mundy
1370 : ! **************************************************************************************************
1371 34545 : SUBROUTINE pw_grid_count(h_inv, pw_grid, mp_comm, cutoff, yz_mask)
1372 :
1373 : REAL(KIND=dp), DIMENSION(3, 3) :: h_inv
1374 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
1375 :
1376 : CLASS(mp_comm_type), INTENT(IN) :: mp_comm
1377 : REAL(KIND=dp), INTENT(IN) :: cutoff
1378 : INTEGER, DIMENSION(:, :), INTENT(OUT) :: yz_mask
1379 :
1380 : INTEGER :: l, m, mm, n, n_upperlimit, nlim(2), nn
1381 : INTEGER(KIND=int_8) :: gpt
1382 : REAL(KIND=dp) :: length, length_x, length_y, length_z
1383 :
1384 : ASSOCIATE (bounds => pw_grid%bounds)
1385 :
1386 34545 : IF (pw_grid%grid_span == HALFSPACE) THEN
1387 : n_upperlimit = 0
1388 31102 : ELSE IF (pw_grid%grid_span == FULLSPACE) THEN
1389 31102 : n_upperlimit = bounds(2, 3)
1390 : ELSE
1391 0 : CPABORT("No type set for the grid")
1392 : END IF
1393 :
1394 : ! finds valid g-points within grid
1395 34545 : gpt = 0
1396 34545 : IF (pw_grid%para%mode == PW_MODE_LOCAL) THEN
1397 5217 : nlim(1) = bounds(1, 3)
1398 5217 : nlim(2) = n_upperlimit
1399 29328 : ELSE IF (pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
1400 29328 : n = n_upperlimit - bounds(1, 3) + 1
1401 29328 : nlim = get_limit(n, mp_comm%num_pe, mp_comm%mepos)
1402 87984 : nlim = nlim + bounds(1, 3) - 1
1403 : ELSE
1404 0 : CPABORT("para % mode not specified")
1405 : END IF
1406 :
1407 33793219 : yz_mask = 0
1408 500204 : DO n = nlim(1), nlim(2)
1409 431114 : nn = n - bounds(1, 3) + 1
1410 16451085 : DO m = bounds(1, 2), bounds(2, 2)
1411 15985426 : mm = m - bounds(1, 2) + 1
1412 873344105 : DO l = bounds(1, 1), bounds(2, 1)
1413 856927565 : IF (pw_grid%grid_span == HALFSPACE .AND. n == 0) THEN
1414 2991708 : IF ((m == 0 .AND. l > 0) .OR. (m > 0)) CYCLE
1415 : END IF
1416 :
1417 855318911 : CALL pw_vec_length(h_inv, length_x, length_y, length_z, length, l, m, n)
1418 :
1419 871304337 : IF (0.5_dp*length <= cutoff) THEN
1420 814382545 : gpt = gpt + 1
1421 814382545 : yz_mask(mm, nn) = yz_mask(mm, nn) + 1
1422 : END IF
1423 :
1424 : END DO
1425 : END DO
1426 : END DO
1427 : END ASSOCIATE
1428 :
1429 : ! number of g-vectors for grid
1430 34545 : IF (pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
1431 29328 : CALL mp_comm%sum(gpt)
1432 62913160 : CALL mp_comm%sum(yz_mask)
1433 : END IF
1434 34545 : pw_grid%ngpts_cut = gpt
1435 :
1436 34545 : END SUBROUTINE pw_grid_count
1437 :
1438 : ! **************************************************************************************************
1439 : !> \brief Setup maps from 1d to 3d space
1440 : !> \param h_inv ...
1441 : !> \param pw_grid ...
1442 : !> \param cutoff ...
1443 : !> \par History
1444 : !> JGH (29-12-2000) : Adapted for parallel use
1445 : !> \author apsi
1446 : !> Christopher Mundy
1447 : ! **************************************************************************************************
1448 34545 : SUBROUTINE pw_grid_assign(h_inv, pw_grid, cutoff)
1449 :
1450 : REAL(KIND=dp), DIMENSION(3, 3) :: h_inv
1451 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
1452 : REAL(KIND=dp), INTENT(IN) :: cutoff
1453 :
1454 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_assign'
1455 :
1456 : INTEGER :: gpt, handle, i, ip, l, lby, lbz, ll, m, &
1457 : mm, n, n_upperlimit, nn
1458 : INTEGER(KIND=int_8) :: gpt_global
1459 : INTEGER, DIMENSION(2, 3) :: bol, bounds
1460 : REAL(KIND=dp) :: length, length_x, length_y, length_z
1461 :
1462 34545 : CALL timeset(routineN, handle)
1463 :
1464 345450 : bounds = pw_grid%bounds
1465 34545 : lby = pw_grid%bounds(1, 2)
1466 34545 : lbz = pw_grid%bounds(1, 3)
1467 :
1468 34545 : IF (pw_grid%grid_span == HALFSPACE) THEN
1469 : n_upperlimit = 0
1470 31102 : ELSE IF (pw_grid%grid_span == FULLSPACE) THEN
1471 31102 : n_upperlimit = bounds(2, 3)
1472 : ELSE
1473 0 : CPABORT("No type set for the grid")
1474 : END IF
1475 :
1476 : ! finds valid g-points within grid
1477 34545 : IF (pw_grid%para%mode == PW_MODE_LOCAL) THEN
1478 5217 : gpt = 0
1479 71713 : DO n = bounds(1, 3), n_upperlimit
1480 1614163 : DO m = bounds(1, 2), bounds(2, 2)
1481 60291842 : DO l = bounds(1, 1), bounds(2, 1)
1482 58682896 : IF (pw_grid%grid_span == HALFSPACE .AND. n == 0) THEN
1483 1162021 : IF ((m == 0 .AND. l > 0) .OR. (m > 0)) CYCLE
1484 : END IF
1485 :
1486 57978184 : CALL pw_vec_length(h_inv, length_x, length_y, length_z, length, l, m, n)
1487 :
1488 59520634 : IF (0.5_dp*length <= cutoff) THEN
1489 44379749 : gpt = gpt + 1
1490 44379749 : pw_grid%g(1, gpt) = length_x
1491 44379749 : pw_grid%g(2, gpt) = length_y
1492 44379749 : pw_grid%g(3, gpt) = length_z
1493 44379749 : pw_grid%gsq(gpt) = length
1494 44379749 : pw_grid%g_hat(1, gpt) = l
1495 44379749 : pw_grid%g_hat(2, gpt) = m
1496 44379749 : pw_grid%g_hat(3, gpt) = n
1497 : END IF
1498 :
1499 : END DO
1500 : END DO
1501 : END DO
1502 :
1503 : ELSE
1504 :
1505 29328 : IF (pw_grid%para%ray_distribution) THEN
1506 :
1507 29316 : gpt = 0
1508 29316 : ip = pw_grid%para%group%mepos
1509 14756568 : DO i = 1, pw_grid%para%nyzray(ip)
1510 14727252 : n = pw_grid%para%yzp(2, i, ip) + lbz - 1
1511 14727252 : m = pw_grid%para%yzp(1, i, ip) + lby - 1
1512 14727252 : IF (n > n_upperlimit) CYCLE
1513 796959329 : DO l = bounds(1, 1), bounds(2, 1)
1514 782805180 : IF (pw_grid%grid_span == HALFSPACE .AND. n == 0) THEN
1515 1649477 : IF ((m == 0 .AND. l > 0) .OR. (m > 0)) CYCLE
1516 : END IF
1517 :
1518 781981329 : CALL pw_vec_length(h_inv, length_x, length_y, length_z, length, l, m, n)
1519 :
1520 796708581 : IF (0.5_dp*length <= cutoff) THEN
1521 769959875 : gpt = gpt + 1
1522 769959875 : pw_grid%g(1, gpt) = length_x
1523 769959875 : pw_grid%g(2, gpt) = length_y
1524 769959875 : pw_grid%g(3, gpt) = length_z
1525 769959875 : pw_grid%gsq(gpt) = length
1526 769959875 : pw_grid%g_hat(1, gpt) = l
1527 769959875 : pw_grid%g_hat(2, gpt) = m
1528 769959875 : pw_grid%g_hat(3, gpt) = n
1529 : END IF
1530 :
1531 : END DO
1532 : END DO
1533 :
1534 : ELSE
1535 :
1536 120 : bol = pw_grid%para%bo(:, :, pw_grid%para%group%mepos, 3)
1537 12 : gpt = 0
1538 254 : DO n = bounds(1, 3), bounds(2, 3)
1539 242 : IF (n < 0) THEN
1540 120 : nn = n + pw_grid%npts(3) + 1
1541 : ELSE
1542 122 : nn = n + 1
1543 : END IF
1544 242 : IF (nn < bol(1, 3) .OR. nn > bol(2, 3)) CYCLE
1545 4816 : DO m = bounds(1, 2), bounds(2, 2)
1546 4562 : IF (m < 0) THEN
1547 2216 : mm = m + pw_grid%npts(2) + 1
1548 : ELSE
1549 2346 : mm = m + 1
1550 : END IF
1551 4562 : IF (mm < bol(1, 2) .OR. mm > bol(2, 2)) CYCLE
1552 45444 : DO l = bounds(1, 1), bounds(2, 1)
1553 42921 : IF (l < 0) THEN
1554 21148 : ll = l + pw_grid%npts(1) + 1
1555 : ELSE
1556 21773 : ll = l + 1
1557 : END IF
1558 42921 : IF (ll < bol(1, 1) .OR. ll > bol(2, 1)) CYCLE
1559 :
1560 42921 : CALL pw_vec_length(h_inv, length_x, length_y, length_z, length, l, m, n)
1561 :
1562 42921 : gpt = gpt + 1
1563 42921 : pw_grid%g(1, gpt) = length_x
1564 42921 : pw_grid%g(2, gpt) = length_y
1565 42921 : pw_grid%g(3, gpt) = length_z
1566 42921 : pw_grid%gsq(gpt) = length
1567 42921 : pw_grid%g_hat(1, gpt) = l
1568 42921 : pw_grid%g_hat(2, gpt) = m
1569 47483 : pw_grid%g_hat(3, gpt) = n
1570 :
1571 : END DO
1572 : END DO
1573 : END DO
1574 :
1575 : END IF
1576 :
1577 : END IF
1578 :
1579 : ! Check the number of g-vectors for grid
1580 34545 : CPASSERT(pw_grid%ngpts_cut_local == gpt)
1581 34545 : IF (pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
1582 29328 : gpt_global = gpt
1583 29328 : CALL pw_grid%para%group%sum(gpt_global)
1584 29328 : CPASSERT(pw_grid%ngpts_cut == gpt_global)
1585 : END IF
1586 :
1587 34545 : pw_grid%have_g0 = .FALSE.
1588 34545 : pw_grid%first_gne0 = 1
1589 613947799 : DO gpt = 1, pw_grid%ngpts_cut_local
1590 625981724 : IF (ALL(pw_grid%g_hat(:, gpt) == 0)) THEN
1591 19881 : pw_grid%have_g0 = .TRUE.
1592 19881 : pw_grid%first_gne0 = 2
1593 19881 : EXIT
1594 : END IF
1595 : END DO
1596 :
1597 : CALL pw_grid_set_maps(pw_grid%grid_span, pw_grid%g_hat, &
1598 34545 : pw_grid%mapl, pw_grid%mapm, pw_grid%mapn, pw_grid%npts)
1599 :
1600 34545 : CALL timestop(handle)
1601 :
1602 34545 : END SUBROUTINE pw_grid_assign
1603 :
1604 : ! **************************************************************************************************
1605 : !> \brief Setup maps from 1d to 3d space
1606 : !> \param grid_span ...
1607 : !> \param g_hat ...
1608 : !> \param mapl ...
1609 : !> \param mapm ...
1610 : !> \param mapn ...
1611 : !> \param npts ...
1612 : !> \par History
1613 : !> JGH (21-12-2000) : Size of g_hat locally determined
1614 : !> \author apsi
1615 : !> Christopher Mundy
1616 : !> \note
1617 : !> Maps are to full 3D space (not distributed)
1618 : ! **************************************************************************************************
1619 34545 : SUBROUTINE pw_grid_set_maps(grid_span, g_hat, mapl, mapm, mapn, npts)
1620 :
1621 : INTEGER, INTENT(IN) :: grid_span
1622 : INTEGER, DIMENSION(:, :), INTENT(IN) :: g_hat
1623 : TYPE(map_pn), INTENT(INOUT) :: mapl, mapm, mapn
1624 : INTEGER, DIMENSION(:), INTENT(IN) :: npts
1625 :
1626 : INTEGER :: gpt, l, m, n, ng
1627 :
1628 34545 : ng = SIZE(g_hat, 2)
1629 :
1630 814417090 : DO gpt = 1, ng
1631 814382545 : l = g_hat(1, gpt)
1632 814382545 : m = g_hat(2, gpt)
1633 814382545 : n = g_hat(3, gpt)
1634 814382545 : IF (l < 0) THEN
1635 404954252 : mapl%pos(l) = l + npts(1)
1636 : ELSE
1637 409428293 : mapl%pos(l) = l
1638 : END IF
1639 814382545 : IF (m < 0) THEN
1640 405373572 : mapm%pos(m) = m + npts(2)
1641 : ELSE
1642 409008973 : mapm%pos(m) = m
1643 : END IF
1644 814382545 : IF (n < 0) THEN
1645 421253873 : mapn%pos(n) = n + npts(3)
1646 : ELSE
1647 393128672 : mapn%pos(n) = n
1648 : END IF
1649 :
1650 : ! Generating the maps to the full 3-d space
1651 :
1652 814417090 : IF (grid_span == HALFSPACE) THEN
1653 :
1654 33501278 : IF (l <= 0) THEN
1655 17243616 : mapl%neg(l) = -l
1656 : ELSE
1657 16257662 : mapl%neg(l) = npts(1) - l
1658 : END IF
1659 33501278 : IF (m <= 0) THEN
1660 17683656 : mapm%neg(m) = -m
1661 : ELSE
1662 15817622 : mapm%neg(m) = npts(2) - m
1663 : END IF
1664 33501278 : IF (n <= 0) THEN
1665 33501278 : mapn%neg(n) = -n
1666 : ELSE
1667 0 : mapn%neg(n) = npts(3) - n
1668 : END IF
1669 :
1670 : END IF
1671 :
1672 : END DO
1673 :
1674 34545 : END SUBROUTINE pw_grid_set_maps
1675 :
1676 : ! **************************************************************************************************
1677 : !> \brief Allocate all (Pointer) Arrays in pw_grid
1678 : !> \param pw_grid ...
1679 : !> \param ng ...
1680 : !> \param bounds ...
1681 : !> \par History
1682 : !> JGH (20-12-2000) : Added status variable
1683 : !> Bounds of arrays now from calling routine, this
1684 : !> makes it independent from parallel setup
1685 : !> \author apsi
1686 : !> Christopher Mundy
1687 : ! **************************************************************************************************
1688 34545 : SUBROUTINE pw_grid_allocate(pw_grid, ng, bounds)
1689 :
1690 : ! Argument
1691 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
1692 : INTEGER, INTENT(IN) :: ng
1693 : INTEGER, DIMENSION(:, :), INTENT(IN) :: bounds
1694 :
1695 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_allocate'
1696 :
1697 : INTEGER :: nmaps
1698 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
1699 : INTEGER(KIND=C_SIZE_T) :: length
1700 : TYPE(C_PTR) :: cptr_g_hatmap
1701 : INTEGER :: stat
1702 : #endif
1703 :
1704 : INTEGER :: handle
1705 :
1706 34545 : CALL timeset(routineN, handle)
1707 :
1708 103635 : ALLOCATE (pw_grid%g(3, ng))
1709 103635 : ALLOCATE (pw_grid%gsq(ng))
1710 103635 : ALLOCATE (pw_grid%g_hat(3, ng))
1711 :
1712 : nmaps = 1
1713 34545 : IF (pw_grid%grid_span == HALFSPACE) nmaps = 2
1714 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
1715 : CALL offload_activate_chosen_device()
1716 :
1717 : length = INT(int_size*MAX(ng, 1)*MAX(nmaps, 1), KIND=C_SIZE_T)
1718 : stat = offload_malloc_pinned_mem(cptr_g_hatmap, length)
1719 : CPASSERT(stat == 0)
1720 : CALL c_f_pointer(cptr_g_hatmap, pw_grid%g_hatmap, (/MAX(ng, 1), MAX(nmaps, 1)/))
1721 : #else
1722 34545 : ALLOCATE (pw_grid%g_hatmap(1, 1))
1723 : #endif
1724 :
1725 34545 : IF (pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
1726 : ALLOCATE (pw_grid%grays(pw_grid%npts(1), &
1727 117312 : pw_grid%para%nyzray(pw_grid%para%group%mepos)))
1728 : END IF
1729 :
1730 103635 : ALLOCATE (pw_grid%mapl%pos(bounds(1, 1):bounds(2, 1)))
1731 103635 : ALLOCATE (pw_grid%mapl%neg(bounds(1, 1):bounds(2, 1)))
1732 103635 : ALLOCATE (pw_grid%mapm%pos(bounds(1, 2):bounds(2, 2)))
1733 103635 : ALLOCATE (pw_grid%mapm%neg(bounds(1, 2):bounds(2, 2)))
1734 103635 : ALLOCATE (pw_grid%mapn%pos(bounds(1, 3):bounds(2, 3)))
1735 103635 : ALLOCATE (pw_grid%mapn%neg(bounds(1, 3):bounds(2, 3)))
1736 :
1737 34545 : CALL timestop(handle)
1738 :
1739 34545 : END SUBROUTINE pw_grid_allocate
1740 :
1741 : ! **************************************************************************************************
1742 : !> \brief Sort g-vectors according to length
1743 : !> \param pw_grid ...
1744 : !> \param ref_grid ...
1745 : !> \par History
1746 : !> JGH (20-12-2000) : allocate idx, ng = SIZE ( pw_grid % gsq ) the
1747 : !> sorting is local and independent from parallelisation
1748 : !> WARNING: Global ordering depends now on the number
1749 : !> of cpus.
1750 : !> JGH (28-02-2001) : check for ordering against reference grid
1751 : !> JGH (01-05-2001) : sort spherical cutoff grids also within shells
1752 : !> reference grids for non-spherical cutoffs
1753 : !> JGH (20-06-2001) : do not sort non-spherical grids
1754 : !> JGH (19-02-2003) : Order all grids, this makes subgrids also for
1755 : !> non-spherical cutoffs possible
1756 : !> JGH (21-02-2003) : Introduce gather array for general reference grids
1757 : !> \author apsi
1758 : !> Christopher Mundy
1759 : ! **************************************************************************************************
1760 34545 : SUBROUTINE pw_grid_sort(pw_grid, ref_grid)
1761 :
1762 : ! Argument
1763 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
1764 : TYPE(pw_grid_type), INTENT(IN), OPTIONAL :: ref_grid
1765 :
1766 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_sort'
1767 :
1768 : INTEGER :: handle, i, ig, ih, ip, is, it, ng, ngr
1769 : INTEGER, ALLOCATABLE, DIMENSION(:) :: idx
1770 34545 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: int_tmp
1771 : LOGICAL :: g_found
1772 : REAL(KIND=dp) :: gig, gigr
1773 34545 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: real_tmp
1774 :
1775 34545 : CALL timeset(routineN, handle)
1776 :
1777 34545 : ng = SIZE(pw_grid%gsq)
1778 103635 : ALLOCATE (idx(ng))
1779 :
1780 : ! grids are (locally) ordered by length of G-vectors
1781 34545 : CALL sort(pw_grid%gsq, ng, idx)
1782 : ! within shells order wrt x,y,z
1783 34545 : CALL sort_shells(pw_grid%gsq, pw_grid%g_hat, idx)
1784 :
1785 103635 : ALLOCATE (real_tmp(3, ng))
1786 814417090 : DO i = 1, ng
1787 814382545 : real_tmp(1, i) = pw_grid%g(1, idx(i))
1788 814382545 : real_tmp(2, i) = pw_grid%g(2, idx(i))
1789 814417090 : real_tmp(3, i) = pw_grid%g(3, idx(i))
1790 : END DO
1791 814417090 : DO i = 1, ng
1792 814382545 : pw_grid%g(1, i) = real_tmp(1, i)
1793 814382545 : pw_grid%g(2, i) = real_tmp(2, i)
1794 814417090 : pw_grid%g(3, i) = real_tmp(3, i)
1795 : END DO
1796 34545 : DEALLOCATE (real_tmp)
1797 :
1798 103635 : ALLOCATE (int_tmp(3, ng))
1799 814417090 : DO i = 1, ng
1800 814382545 : int_tmp(1, i) = pw_grid%g_hat(1, idx(i))
1801 814382545 : int_tmp(2, i) = pw_grid%g_hat(2, idx(i))
1802 814417090 : int_tmp(3, i) = pw_grid%g_hat(3, idx(i))
1803 : END DO
1804 814417090 : DO i = 1, ng
1805 814382545 : pw_grid%g_hat(1, i) = int_tmp(1, i)
1806 814382545 : pw_grid%g_hat(2, i) = int_tmp(2, i)
1807 814417090 : pw_grid%g_hat(3, i) = int_tmp(3, i)
1808 : END DO
1809 34545 : DEALLOCATE (int_tmp)
1810 :
1811 34545 : DEALLOCATE (idx)
1812 :
1813 : ! check if ordering is compatible to reference grid
1814 34545 : IF (PRESENT(ref_grid)) THEN
1815 21076 : ngr = SIZE(ref_grid%gsq)
1816 21076 : ngr = MIN(ng, ngr)
1817 21076 : IF (pw_grid%spherical) THEN
1818 0 : IF (.NOT. ALL(pw_grid%g_hat(1:3, 1:ngr) &
1819 : == ref_grid%g_hat(1:3, 1:ngr))) THEN
1820 0 : CPABORT("G space sorting not compatible")
1821 : END IF
1822 : ELSE
1823 63228 : ALLOCATE (pw_grid%gidx(1:ngr))
1824 179122307 : pw_grid%gidx = 0
1825 : ! first try as many trivial associations as possible
1826 21076 : it = 0
1827 95137446 : DO ig = 1, ngr
1828 380489381 : IF (.NOT. ALL(pw_grid%g_hat(1:3, ig) &
1829 : == ref_grid%g_hat(1:3, ig))) EXIT
1830 95116370 : pw_grid%gidx(ig) = ig
1831 95137446 : it = ig
1832 : END DO
1833 : ! now for the ones that could not be done
1834 21076 : IF (ng == ngr) THEN
1835 : ! for the case pw_grid <= ref_grid
1836 21076 : is = it
1837 84005937 : DO ig = it + 1, ngr
1838 83984861 : gig = pw_grid%gsq(ig)
1839 83984861 : gigr = MAX(1._dp, gig)
1840 83984861 : g_found = .FALSE.
1841 420634785 : DO ih = is + 1, SIZE(ref_grid%gsq)
1842 420634785 : IF (ABS(gig - ref_grid%gsq(ih))/gigr > 1.e-12_dp) CYCLE
1843 : g_found = .TRUE.
1844 336649924 : EXIT
1845 : END DO
1846 83984861 : IF (.NOT. g_found) THEN
1847 0 : WRITE (*, "(A,I10,F20.10)") "G-vector", ig, pw_grid%gsq(ig)
1848 0 : CPABORT("G vector not found")
1849 : END IF
1850 83984861 : ip = ih - 1
1851 5797599866 : DO ih = ip + 1, SIZE(ref_grid%gsq)
1852 5797599866 : IF (ABS(gig - ref_grid%gsq(ih))/gigr > 1.e-12_dp) CYCLE
1853 5797599866 : IF (pw_grid%g_hat(1, ig) /= ref_grid%g_hat(1, ih)) CYCLE
1854 320676259 : IF (pw_grid%g_hat(2, ig) /= ref_grid%g_hat(2, ih)) CYCLE
1855 103942135 : IF (pw_grid%g_hat(3, ig) /= ref_grid%g_hat(3, ih)) CYCLE
1856 83984861 : pw_grid%gidx(ig) = ih
1857 5797599866 : EXIT
1858 : END DO
1859 83984861 : IF (pw_grid%gidx(ig) == 0) THEN
1860 0 : WRITE (*, "(A,2I10)") " G-Shell ", is + 1, ip + 1
1861 : WRITE (*, "(A,I10,3I6,F20.10)") &
1862 0 : " G-vector", ig, pw_grid%g_hat(1:3, ig), pw_grid%gsq(ig)
1863 0 : DO ih = 1, SIZE(ref_grid%gsq)
1864 0 : IF (pw_grid%g_hat(1, ig) /= ref_grid%g_hat(1, ih)) CYCLE
1865 0 : IF (pw_grid%g_hat(2, ig) /= ref_grid%g_hat(2, ih)) CYCLE
1866 0 : IF (pw_grid%g_hat(3, ig) /= ref_grid%g_hat(3, ih)) CYCLE
1867 : WRITE (*, "(A,I10,3I6,F20.10)") &
1868 0 : " G-vector", ih, ref_grid%g_hat(1:3, ih), ref_grid%gsq(ih)
1869 : END DO
1870 0 : CPABORT("G vector not found")
1871 : END IF
1872 84005937 : is = ip
1873 : END DO
1874 : ELSE
1875 : ! for the case pw_grid > ref_grid
1876 0 : is = it
1877 0 : DO ig = it + 1, ngr
1878 0 : gig = ref_grid%gsq(ig)
1879 0 : gigr = MAX(1._dp, gig)
1880 0 : g_found = .FALSE.
1881 0 : DO ih = is + 1, ng
1882 0 : IF (ABS(pw_grid%gsq(ih) - gig)/gigr > 1.e-12_dp) CYCLE
1883 : g_found = .TRUE.
1884 0 : EXIT
1885 : END DO
1886 0 : IF (.NOT. g_found) THEN
1887 0 : WRITE (*, "(A,I10,F20.10)") "G-vector", ig, ref_grid%gsq(ig)
1888 0 : CPABORT("G vector not found")
1889 : END IF
1890 0 : ip = ih - 1
1891 0 : DO ih = ip + 1, ng
1892 0 : IF (ABS(pw_grid%gsq(ih) - gig)/gigr > 1.e-12_dp) CYCLE
1893 0 : IF (pw_grid%g_hat(1, ih) /= ref_grid%g_hat(1, ig)) CYCLE
1894 0 : IF (pw_grid%g_hat(2, ih) /= ref_grid%g_hat(2, ig)) CYCLE
1895 0 : IF (pw_grid%g_hat(3, ih) /= ref_grid%g_hat(3, ig)) CYCLE
1896 0 : pw_grid%gidx(ig) = ih
1897 0 : EXIT
1898 : END DO
1899 0 : IF (pw_grid%gidx(ig) == 0) THEN
1900 0 : WRITE (*, "(A,2I10)") " G-Shell ", is + 1, ip + 1
1901 : WRITE (*, "(A,I10,3I6,F20.10)") &
1902 0 : " G-vector", ig, ref_grid%g_hat(1:3, ig), ref_grid%gsq(ig)
1903 0 : DO ih = 1, ng
1904 0 : IF (pw_grid%g_hat(1, ih) /= ref_grid%g_hat(1, ig)) CYCLE
1905 0 : IF (pw_grid%g_hat(2, ih) /= ref_grid%g_hat(2, ig)) CYCLE
1906 0 : IF (pw_grid%g_hat(3, ih) /= ref_grid%g_hat(3, ig)) CYCLE
1907 : WRITE (*, "(A,I10,3I6,F20.10)") &
1908 0 : " G-vector", ih, pw_grid%g_hat(1:3, ih), pw_grid%gsq(ih)
1909 : END DO
1910 0 : CPABORT("G vector not found")
1911 : END IF
1912 0 : is = ip
1913 : END DO
1914 : END IF
1915 : ! test if all g-vectors are associated
1916 179122307 : IF (ANY(pw_grid%gidx == 0)) THEN
1917 0 : CPABORT("G space sorting not compatible")
1918 : END IF
1919 : END IF
1920 : END IF
1921 :
1922 : !check if G=0 is at first position
1923 34545 : IF (pw_grid%have_g0) THEN
1924 19881 : IF (pw_grid%gsq(1) /= 0.0_dp) THEN
1925 0 : CPABORT("G = 0 not in first position")
1926 : END IF
1927 : END IF
1928 :
1929 34545 : CALL timestop(handle)
1930 :
1931 34545 : END SUBROUTINE pw_grid_sort
1932 :
1933 : ! **************************************************************************************************
1934 : !> \brief ...
1935 : !> \param gsq ...
1936 : !> \param g_hat ...
1937 : !> \param idx ...
1938 : ! **************************************************************************************************
1939 34545 : SUBROUTINE sort_shells(gsq, g_hat, idx)
1940 :
1941 : ! Argument
1942 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: gsq
1943 : INTEGER, DIMENSION(:, :), INTENT(IN) :: g_hat
1944 : INTEGER, DIMENSION(:), INTENT(INOUT) :: idx
1945 :
1946 : CHARACTER(len=*), PARAMETER :: routineN = 'sort_shells'
1947 : REAL(KIND=dp), PARAMETER :: small = 5.e-16_dp
1948 :
1949 : INTEGER :: handle, ig, ng, s1, s2
1950 : REAL(KIND=dp) :: s_begin
1951 :
1952 34545 : CALL timeset(routineN, handle)
1953 :
1954 : ! Juergs temporary hack to get the grids sorted for large (4000Ry) cutoffs.
1955 : ! might need to call lapack for machine precision.
1956 :
1957 34545 : ng = SIZE(gsq)
1958 34545 : s_begin = -1.0_dp
1959 34545 : s1 = 0
1960 34545 : s2 = 0
1961 : ig = 1
1962 814417090 : DO ig = 1, ng
1963 814417090 : IF (ABS(gsq(ig) - s_begin) < small) THEN
1964 772760147 : s2 = ig
1965 : ELSE
1966 41622398 : CALL redist(g_hat, idx, s1, s2)
1967 41622398 : s_begin = gsq(ig)
1968 41622398 : s1 = ig
1969 41622398 : s2 = ig
1970 : END IF
1971 : END DO
1972 34545 : CALL redist(g_hat, idx, s1, s2)
1973 :
1974 34545 : CALL timestop(handle)
1975 :
1976 34545 : END SUBROUTINE sort_shells
1977 :
1978 : ! **************************************************************************************************
1979 : !> \brief ...
1980 : !> \param g_hat ...
1981 : !> \param idx ...
1982 : !> \param s1 ...
1983 : !> \param s2 ...
1984 : ! **************************************************************************************************
1985 41656943 : SUBROUTINE redist(g_hat, idx, s1, s2)
1986 :
1987 : ! Argument
1988 : INTEGER, DIMENSION(:, :), INTENT(IN) :: g_hat
1989 : INTEGER, DIMENSION(:), INTENT(INOUT) :: idx
1990 : INTEGER, INTENT(IN) :: s1, s2
1991 :
1992 : INTEGER :: i, ii, n1, n2, n3, ns
1993 41656943 : INTEGER, ALLOCATABLE, DIMENSION(:) :: indl
1994 41656943 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: slen
1995 :
1996 41656943 : IF (s2 <= s1) RETURN
1997 39774203 : ns = s2 - s1 + 1
1998 119322609 : ALLOCATE (indl(ns))
1999 119322609 : ALLOCATE (slen(ns))
2000 :
2001 852308553 : DO i = s1, s2
2002 812534350 : ii = idx(i)
2003 812534350 : n1 = g_hat(1, ii)
2004 812534350 : n2 = g_hat(2, ii)
2005 812534350 : n3 = g_hat(3, ii)
2006 : slen(i - s1 + 1) = 1000.0_dp*REAL(n1, dp) + &
2007 852308553 : REAL(n2, dp) + 0.001_dp*REAL(n3, dp)
2008 : END DO
2009 39774203 : CALL sort(slen, ns, indl)
2010 852308553 : DO i = 1, ns
2011 812534350 : ii = indl(i) + s1 - 1
2012 852308553 : indl(i) = idx(ii)
2013 : END DO
2014 852308553 : idx(s1:s2) = indl(1:ns)
2015 :
2016 39774203 : DEALLOCATE (indl)
2017 39774203 : DEALLOCATE (slen)
2018 :
2019 : END SUBROUTINE redist
2020 :
2021 : ! **************************************************************************************************
2022 : !> \brief Reorder yzq and yzp arrays for parallel FFT according to FFT mapping
2023 : !> \param pw_grid ...
2024 : !> \param yz ...
2025 : !> \par History
2026 : !> none
2027 : !> \author JGH (17-Jan-2001)
2028 : ! **************************************************************************************************
2029 34545 : SUBROUTINE pw_grid_remap(pw_grid, yz)
2030 :
2031 : ! Argument
2032 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
2033 : INTEGER, DIMENSION(:, :), INTENT(OUT) :: yz
2034 :
2035 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_remap'
2036 :
2037 : INTEGER :: gpt, handle, i, ip, is, j, m, n
2038 :
2039 34545 : IF (pw_grid%para%mode == PW_MODE_LOCAL) RETURN
2040 :
2041 29328 : CALL timeset(routineN, handle)
2042 :
2043 : ASSOCIATE (ny => pw_grid%npts(2), nz => pw_grid%npts(3), posm => pw_grid%mapm%pos, posn => pw_grid%mapn%pos, &
2044 : negm => pw_grid%mapm%neg, negn => pw_grid%mapn%neg)
2045 :
2046 31471244 : yz = 0
2047 88562244 : pw_grid%para%yzp = 0
2048 31471244 : pw_grid%para%yzq = 0
2049 :
2050 770032124 : DO gpt = 1, SIZE(pw_grid%gsq)
2051 770002796 : m = posm(pw_grid%g_hat(2, gpt)) + 1
2052 770002796 : n = posn(pw_grid%g_hat(3, gpt)) + 1
2053 770032124 : yz(m, n) = yz(m, n) + 1
2054 : END DO
2055 29328 : IF (pw_grid%grid_span == HALFSPACE) THEN
2056 20940275 : DO gpt = 1, SIZE(pw_grid%gsq)
2057 20937985 : m = negm(pw_grid%g_hat(2, gpt)) + 1
2058 20937985 : n = negn(pw_grid%g_hat(3, gpt)) + 1
2059 20940275 : yz(m, n) = yz(m, n) + 1
2060 : END DO
2061 : END IF
2062 :
2063 29328 : ip = pw_grid%para%group%mepos
2064 29328 : is = 0
2065 827966 : DO i = 1, nz
2066 31471244 : DO j = 1, ny
2067 31441916 : IF (yz(j, i) > 0) THEN
2068 14729533 : is = is + 1
2069 14729533 : pw_grid%para%yzp(1, is, ip) = j
2070 14729533 : pw_grid%para%yzp(2, is, ip) = i
2071 14729533 : pw_grid%para%yzq(j, i) = is
2072 : END IF
2073 : END DO
2074 : END DO
2075 : END ASSOCIATE
2076 :
2077 29328 : CPASSERT(is == pw_grid%para%nyzray(ip))
2078 29328 : CALL pw_grid%para%group%sum(pw_grid%para%yzp)
2079 :
2080 29328 : CALL timestop(handle)
2081 :
2082 29328 : END SUBROUTINE pw_grid_remap
2083 :
2084 : ! **************************************************************************************************
2085 : !> \brief Recalculate the g-vectors after a change of the box
2086 : !> \param cell_hmat ...
2087 : !> \param pw_grid ...
2088 : !> \par History
2089 : !> JGH (20-12-2000) : get local grid size from definition of g.
2090 : !> Assume that gsq is allocated.
2091 : !> Local routine, no information on distribution of
2092 : !> PW required.
2093 : !> JGH (8-Mar-2001) : also update information on volume and grid spaceing
2094 : !> \author apsi
2095 : !> Christopher Mundy
2096 : ! **************************************************************************************************
2097 11416 : SUBROUTINE pw_grid_change(cell_hmat, pw_grid)
2098 : ! Argument
2099 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: cell_hmat
2100 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
2101 :
2102 : INTEGER :: gpt
2103 : REAL(KIND=dp) :: cell_deth, l, m, n
2104 : REAL(KIND=dp), DIMENSION(3, 3) :: cell_h_inv
2105 11416 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: g
2106 :
2107 11416 : cell_deth = ABS(det_3x3(cell_hmat))
2108 11416 : IF (cell_deth < 1.0E-10_dp) THEN
2109 : CALL cp_abort(__LOCATION__, &
2110 : "An invalid set of cell vectors was specified. "// &
2111 0 : "The determinant det(h) is too small")
2112 : END IF
2113 11416 : cell_h_inv = inv_3x3(cell_hmat)
2114 :
2115 11416 : g => pw_grid%g
2116 :
2117 11416 : CALL cell2grid(cell_hmat, cell_h_inv, cell_deth, pw_grid)
2118 :
2119 105739451 : DO gpt = 1, SIZE(g, 2)
2120 :
2121 105728035 : l = twopi*REAL(pw_grid%g_hat(1, gpt), KIND=dp)
2122 105728035 : m = twopi*REAL(pw_grid%g_hat(2, gpt), KIND=dp)
2123 105728035 : n = twopi*REAL(pw_grid%g_hat(3, gpt), KIND=dp)
2124 :
2125 105728035 : g(1, gpt) = l*cell_h_inv(1, 1) + m*cell_h_inv(2, 1) + n*cell_h_inv(3, 1)
2126 105728035 : g(2, gpt) = l*cell_h_inv(1, 2) + m*cell_h_inv(2, 2) + n*cell_h_inv(3, 2)
2127 105728035 : g(3, gpt) = l*cell_h_inv(1, 3) + m*cell_h_inv(2, 3) + n*cell_h_inv(3, 3)
2128 :
2129 : pw_grid%gsq(gpt) = g(1, gpt)*g(1, gpt) &
2130 : + g(2, gpt)*g(2, gpt) &
2131 105739451 : + g(3, gpt)*g(3, gpt)
2132 :
2133 : END DO
2134 :
2135 11416 : END SUBROUTINE pw_grid_change
2136 :
2137 : ! **************************************************************************************************
2138 : !> \brief retains the given pw grid
2139 : !> \param pw_grid the pw grid to retain
2140 : !> \par History
2141 : !> 04.2003 created [fawzi]
2142 : !> \author fawzi
2143 : !> \note
2144 : !> see doc/ReferenceCounting.html
2145 : ! **************************************************************************************************
2146 7435678 : SUBROUTINE pw_grid_retain(pw_grid)
2147 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
2148 :
2149 7435678 : CPASSERT(pw_grid%ref_count > 0)
2150 7435678 : pw_grid%ref_count = pw_grid%ref_count + 1
2151 7435678 : END SUBROUTINE pw_grid_retain
2152 :
2153 : ! **************************************************************************************************
2154 : !> \brief releases the given pw grid
2155 : !> \param pw_grid the pw grid to release
2156 : !> \par History
2157 : !> 04.2003 created [fawzi]
2158 : !> \author fawzi
2159 : !> \note
2160 : !> see doc/ReferenceCounting.html
2161 : ! **************************************************************************************************
2162 15298197 : SUBROUTINE pw_grid_release(pw_grid)
2163 :
2164 : TYPE(pw_grid_type), POINTER :: pw_grid
2165 :
2166 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2167 : INTEGER, POINTER :: dummy_ptr
2168 : INTEGER :: stat
2169 : #endif
2170 :
2171 15298197 : IF (ASSOCIATED(pw_grid)) THEN
2172 7531499 : CPASSERT(pw_grid%ref_count > 0)
2173 7531499 : pw_grid%ref_count = pw_grid%ref_count - 1
2174 7531499 : IF (pw_grid%ref_count == 0) THEN
2175 95821 : IF (ASSOCIATED(pw_grid%gidx)) THEN
2176 21076 : DEALLOCATE (pw_grid%gidx)
2177 : END IF
2178 95821 : IF (ASSOCIATED(pw_grid%g)) THEN
2179 34545 : DEALLOCATE (pw_grid%g)
2180 : END IF
2181 95821 : IF (ASSOCIATED(pw_grid%gsq)) THEN
2182 34545 : DEALLOCATE (pw_grid%gsq)
2183 : END IF
2184 95821 : IF (ALLOCATED(pw_grid%g_hat)) THEN
2185 34545 : DEALLOCATE (pw_grid%g_hat)
2186 : END IF
2187 95821 : IF (ASSOCIATED(pw_grid%g_hatmap)) THEN
2188 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2189 : dummy_ptr => pw_grid%g_hatmap(1, 1)
2190 : stat = offload_free_pinned_mem(c_loc(dummy_ptr))
2191 : CPASSERT(stat == 0)
2192 : #else
2193 34545 : DEALLOCATE (pw_grid%g_hatmap)
2194 : #endif
2195 : END IF
2196 95821 : IF (ASSOCIATED(pw_grid%grays)) THEN
2197 29328 : DEALLOCATE (pw_grid%grays)
2198 : END IF
2199 95821 : IF (ALLOCATED(pw_grid%mapl%pos)) THEN
2200 34545 : DEALLOCATE (pw_grid%mapl%pos)
2201 : END IF
2202 95821 : IF (ALLOCATED(pw_grid%mapm%pos)) THEN
2203 34545 : DEALLOCATE (pw_grid%mapm%pos)
2204 : END IF
2205 95821 : IF (ALLOCATED(pw_grid%mapn%pos)) THEN
2206 34545 : DEALLOCATE (pw_grid%mapn%pos)
2207 : END IF
2208 95821 : IF (ALLOCATED(pw_grid%mapl%neg)) THEN
2209 34545 : DEALLOCATE (pw_grid%mapl%neg)
2210 : END IF
2211 95821 : IF (ALLOCATED(pw_grid%mapm%neg)) THEN
2212 34545 : DEALLOCATE (pw_grid%mapm%neg)
2213 : END IF
2214 95821 : IF (ALLOCATED(pw_grid%mapn%neg)) THEN
2215 34545 : DEALLOCATE (pw_grid%mapn%neg)
2216 : END IF
2217 95821 : IF (ALLOCATED(pw_grid%para%bo)) THEN
2218 34545 : DEALLOCATE (pw_grid%para%bo)
2219 : END IF
2220 95821 : IF (pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
2221 29336 : IF (ALLOCATED(pw_grid%para%yzp)) THEN
2222 29328 : DEALLOCATE (pw_grid%para%yzp)
2223 : END IF
2224 29336 : IF (ALLOCATED(pw_grid%para%yzq)) THEN
2225 29328 : DEALLOCATE (pw_grid%para%yzq)
2226 : END IF
2227 29336 : IF (ALLOCATED(pw_grid%para%nyzray)) THEN
2228 29328 : DEALLOCATE (pw_grid%para%nyzray)
2229 : END IF
2230 : END IF
2231 : ! also release groups
2232 95821 : CALL pw_grid%para%group%free()
2233 95821 : IF (ALLOCATED(pw_grid%para%pos_of_x)) THEN
2234 35553 : DEALLOCATE (pw_grid%para%pos_of_x)
2235 : END IF
2236 :
2237 95821 : IF (ASSOCIATED(pw_grid)) THEN
2238 95821 : DEALLOCATE (pw_grid)
2239 : END IF
2240 : END IF
2241 : END IF
2242 15298197 : NULLIFY (pw_grid)
2243 15298197 : END SUBROUTINE pw_grid_release
2244 :
2245 : END MODULE pw_grids
|
