Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 59066 : 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 59066 : CPASSERT(.NOT. ASSOCIATED(pw_grid))
104 3484894 : ALLOCATE (pw_grid)
105 590660 : pw_grid%bounds = bounds
106 590660 : pw_grid%bounds_local = bounds
107 236264 : pw_grid%npts = bounds(2, :) - bounds(1, :) + 1
108 236264 : pw_grid%npts_local = pw_grid%npts
109 236264 : pw_grid%ngpts = PRODUCT(INT(pw_grid%npts, KIND=int_8))
110 59066 : pw_grid%ngpts_cut = pw_grid%ngpts
111 236264 : pw_grid%ngpts_local = PRODUCT(pw_grid%npts)
112 59066 : pw_grid%ngpts_cut_local = pw_grid%ngpts_local
113 59066 : pw_grid%grid_span = FULLSPACE
114 59066 : pw_grid%para%mode = PW_MODE_LOCAL
115 59066 : pw_grid%reference = 0
116 59066 : pw_grid%ref_count = 1
117 59066 : NULLIFY (pw_grid%g)
118 59066 : NULLIFY (pw_grid%gsq)
119 59066 : NULLIFY (pw_grid%g_hatmap)
120 59066 : NULLIFY (pw_grid%gidx)
121 59066 : NULLIFY (pw_grid%grays)
122 :
123 : ! assign a unique tag to this grid
124 59066 : grid_tag = grid_tag + 1
125 59066 : pw_grid%id_nr = grid_tag
126 :
127 : ! parallel info
128 177198 : rs_dims = 1
129 59066 : CALL pw_grid%para%group%create(mp_comm_self, 2, rs_dims)
130 59066 : IF (pw_grid%para%group%num_pe > 1) THEN
131 0 : pw_grid%para%mode = PW_MODE_DISTRIBUTED
132 : ELSE
133 59066 : pw_grid%para%mode = PW_MODE_LOCAL
134 : END IF
135 :
136 59066 : 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 7174059 : 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 7174059 : 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 7174059 : 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 60482 : 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 60482 : CPASSERT(pw_grid%ref_count > 0)
198 :
199 60482 : IF (PRESENT(id_nr)) id_nr = pw_grid%id_nr
200 60482 : IF (PRESENT(mode)) mode = pw_grid%para%mode
201 60482 : IF (PRESENT(vol)) vol = pw_grid%vol
202 60482 : IF (PRESENT(dvol)) dvol = pw_grid%dvol
203 302234 : IF (PRESENT(npts)) npts(1:3) = pw_grid%npts(1:3)
204 60482 : IF (PRESENT(ngpts)) ngpts = pw_grid%ngpts
205 60482 : IF (PRESENT(ngpts_cut)) ngpts_cut = pw_grid%ngpts_cut
206 60618 : IF (PRESENT(dr)) dr = pw_grid%dr
207 60482 : IF (PRESENT(cutoff)) cutoff = pw_grid%cutoff
208 60482 : IF (PRESENT(orthorhombic)) orthorhombic = pw_grid%orthorhombic
209 60482 : IF (PRESENT(gvectors)) gvectors => pw_grid%g
210 60482 : IF (PRESENT(gsquare)) gsquare => pw_grid%gsq
211 :
212 60482 : 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 65160 : 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 65160 : CPASSERT(pw_grid%ref_count > 0)
236 :
237 65160 : IF (PRESENT(grid_span)) THEN
238 32013 : pw_grid%grid_span = grid_span
239 : END IF
240 65160 : 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 65160 : ELSE IF (PRESENT(bounds)) THEN
245 11170 : pw_grid%bounds = bounds
246 4468 : pw_grid%npts = bounds(2, :) - bounds(1, :) + 1
247 64043 : ELSE IF (PRESENT(npts)) THEN
248 128120 : pw_grid%npts = npts
249 320300 : pw_grid%bounds = pw_grid_bounds_from_n(npts)
250 : END IF
251 65160 : IF (PRESENT(cutoff)) THEN
252 33147 : pw_grid%cutoff = cutoff
253 33147 : IF (PRESENT(spherical)) THEN
254 33147 : pw_grid%spherical = spherical
255 : ELSE
256 0 : pw_grid%spherical = .FALSE.
257 : END IF
258 : END IF
259 :
260 65160 : 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 33147 : 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 33147 : CALL timeset(routineN, handle)
313 :
314 33147 : CPASSERT(.NOT. ASSOCIATED(pw_grid))
315 1723644 : ALLOCATE (pw_grid)
316 331470 : pw_grid%bounds = 0
317 33147 : pw_grid%cutoff = 0.0_dp
318 33147 : pw_grid%grid_span = FULLSPACE
319 : pw_grid%para%mode = PW_MODE_LOCAL
320 : pw_grid%reference = 0
321 33147 : pw_grid%ref_count = 1
322 33147 : NULLIFY (pw_grid%g)
323 33147 : NULLIFY (pw_grid%gsq)
324 33147 : NULLIFY (pw_grid%g_hatmap)
325 33147 : NULLIFY (pw_grid%gidx)
326 33147 : NULLIFY (pw_grid%grays)
327 :
328 : ! assign a unique tag to this grid
329 33147 : grid_tag = grid_tag + 1
330 33147 : pw_grid%id_nr = grid_tag
331 :
332 : ! parallel info
333 33147 : IF (mp_comm%num_pe > 1) THEN
334 28044 : pw_grid%para%mode = PW_MODE_DISTRIBUTED
335 : ELSE
336 : pw_grid%para%mode = PW_MODE_LOCAL
337 : END IF
338 :
339 33147 : cell_deth = ABS(det_3x3(cell_hmat))
340 33147 : 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 33147 : cell_h_inv = inv_3x3(cell_hmat)
346 :
347 33147 : IF (PRESENT(grid_span)) THEN
348 32013 : CALL set_pw_grid_info(pw_grid, grid_span=grid_span)
349 : END IF
350 :
351 33147 : IF (PRESENT(spherical)) THEN
352 32935 : my_spherical = spherical
353 : ELSE
354 212 : my_spherical = .FALSE.
355 : END IF
356 :
357 33147 : IF (PRESENT(odd)) THEN
358 28632 : my_odd = odd
359 : ELSE
360 4515 : my_odd = .FALSE.
361 : END IF
362 :
363 33147 : IF (PRESENT(fft_usage)) THEN
364 33023 : my_fft_usage = fft_usage
365 : ELSE
366 124 : my_fft_usage = .FALSE.
367 : END IF
368 :
369 33147 : IF (PRESENT(ncommensurate)) THEN
370 28572 : my_ncommensurate = ncommensurate
371 28572 : IF (PRESENT(icommensurate)) THEN
372 28572 : my_icommensurate = icommensurate
373 : ELSE
374 0 : my_icommensurate = MIN(1, ncommensurate)
375 : END IF
376 : ELSE
377 4575 : my_ncommensurate = 0
378 4575 : my_icommensurate = 1
379 : END IF
380 :
381 33147 : IF (PRESENT(bounds)) THEN
382 1117 : IF (PRESENT(cutoff)) THEN
383 : CALL set_pw_grid_info(pw_grid, bounds=bounds, cutoff=cutoff, &
384 172 : spherical=my_spherical)
385 : ELSE
386 3780 : n = bounds(2, :) - bounds(1, :) + 1
387 945 : my_cutoff = pw_find_cutoff(n, cell_h_inv)
388 945 : my_cutoff = 0.5_dp*my_cutoff*my_cutoff
389 : CALL set_pw_grid_info(pw_grid, bounds=bounds, cutoff=my_cutoff, &
390 945 : spherical=my_spherical)
391 : END IF
392 32030 : ELSE IF (PRESENT(npts)) THEN
393 2368 : n = npts
394 2368 : IF (PRESENT(cutoff)) THEN
395 22 : my_cutoff = cutoff
396 : ELSE
397 2346 : my_cutoff = pw_find_cutoff(npts, cell_h_inv)
398 2346 : my_cutoff = 0.5_dp*my_cutoff*my_cutoff
399 : END IF
400 2368 : 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 9472 : ref_grid=ref_grid, n_orig=n)
405 : END IF
406 : CALL set_pw_grid_info(pw_grid, npts=n, cutoff=my_cutoff, &
407 2368 : spherical=my_spherical)
408 29662 : 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 29662 : ref_grid=ref_grid)
413 : CALL set_pw_grid_info(pw_grid, npts=n, cutoff=cutoff, &
414 29662 : 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 33147 : 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 33147 : CALL timestop(handle)
427 :
428 33147 : 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 33147 : 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 33147 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: yz_mask
550 : REAL(KIND=dp) :: ecut
551 :
552 : !------------------------------------------------------------------------------
553 :
554 33147 : CALL timeset(routineN, handle)
555 :
556 33147 : CPASSERT(pw_grid%ref_count > 0)
557 :
558 : ! set pointer to possible reference grid
559 33147 : IF (PRESENT(ref_grid)) THEN
560 20912 : pw_grid%reference = ref_grid%id_nr
561 : END IF
562 :
563 33147 : IF (pw_grid%spherical) THEN
564 3207 : ecut = pw_grid%cutoff
565 : ELSE
566 29940 : ecut = 1.e10_dp
567 : END IF
568 :
569 132588 : 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 132588 : ALLOCATE (yz_mask(n(2), n(3)))
576 33147 : CALL pw_grid_count(cell_h_inv, pw_grid, mp_comm, ecut, yz_mask)
577 :
578 : ! Check if reference grid is compatible
579 33147 : IF (PRESENT(ref_grid)) THEN
580 20912 : CPASSERT(pw_grid%para%mode == ref_grid%para%mode)
581 20912 : CPASSERT(pw_grid%grid_span == ref_grid%grid_span)
582 20912 : 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 33147 : rs_dims=rs_dims)
588 :
589 : ! Allocate the grid fields
590 : CALL pw_grid_allocate(pw_grid, pw_grid%ngpts_cut_local, &
591 33147 : pw_grid%bounds)
592 :
593 : ! Fill in the grid structure
594 33147 : CALL pw_grid_assign(cell_h_inv, pw_grid, ecut)
595 :
596 : ! Sort g vector wrt length (only local for each processor)
597 33147 : CALL pw_grid_sort(pw_grid, ref_grid)
598 :
599 33147 : CALL pw_grid_remap(pw_grid, yz_mask)
600 :
601 33147 : DEALLOCATE (yz_mask)
602 :
603 33147 : CALL cell2grid(cell_hmat, cell_h_inv, cell_deth, pw_grid)
604 : !
605 : ! Output: All the information of this grid type
606 : !
607 :
608 33147 : IF (PRESENT(iounit)) THEN
609 33125 : CALL pw_grid_print(pw_grid, iounit)
610 : END IF
611 :
612 33147 : CALL timestop(handle)
613 :
614 33147 : 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 44261 : 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 44261 : pw_grid%vol = ABS(cell_deth)
631 44261 : pw_grid%dvol = pw_grid%vol/REAL(pw_grid%ngpts, KIND=dp)
632 : pw_grid%dr(1) = SQRT(SUM(cell_hmat(:, 1)**2)) &
633 177044 : /REAL(pw_grid%npts(1), KIND=dp)
634 : pw_grid%dr(2) = SQRT(SUM(cell_hmat(:, 2)**2)) &
635 177044 : /REAL(pw_grid%npts(2), KIND=dp)
636 : pw_grid%dr(3) = SQRT(SUM(cell_hmat(:, 3)**2)) &
637 177044 : /REAL(pw_grid%npts(3), KIND=dp)
638 177044 : pw_grid%dh(:, 1) = cell_hmat(:, 1)/REAL(pw_grid%npts(1), KIND=dp)
639 177044 : pw_grid%dh(:, 2) = cell_hmat(:, 2)/REAL(pw_grid%npts(2), KIND=dp)
640 177044 : pw_grid%dh(:, 3) = cell_hmat(:, 3)/REAL(pw_grid%npts(3), KIND=dp)
641 177044 : pw_grid%dh_inv(1, :) = cell_h_inv(1, :)*REAL(pw_grid%npts(1), KIND=dp)
642 177044 : pw_grid%dh_inv(2, :) = cell_h_inv(2, :)*REAL(pw_grid%npts(2), KIND=dp)
643 177044 : 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 44261 : (cell_hmat(3, 1) == 0.0_dp) .AND. (cell_hmat(3, 2) == 0.0_dp)) THEN
648 36555 : pw_grid%orthorhombic = .TRUE.
649 : ELSE
650 7706 : pw_grid%orthorhombic = .FALSE.
651 : END IF
652 44261 : 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 33125 : 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 33125 : IF (pw_grid%para%mode == PW_MODE_LOCAL) THEN
675 5103 : IF (info > 0) THEN
676 : WRITE (info, '(/,A,T71,I10)') &
677 429 : " PW_GRID| Information for grid number ", pw_grid%id_nr
678 429 : 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 429 : WRITE (info, '(" PW_GRID| Cutoff [a.u.]",T71,f10.1)') pw_grid%cutoff
683 429 : 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 229 : WRITE (info, '(A,T78,A)') " PW_GRID| spherical cutoff: ", " NO"
689 : END IF
690 1716 : DO i = 1, 3
691 1287 : WRITE (info, '(A,I3,T30,2I8,T62,A,T71,I10)') " PW_GRID| Bounds ", &
692 1287 : i, pw_grid%bounds(1, I), pw_grid%bounds(2, I), &
693 3003 : "Points:", pw_grid%npts(I)
694 : END DO
695 : WRITE (info, '(A,G12.4,T50,A,T67,F14.4)') &
696 429 : " PW_GRID| Volume element (a.u.^3)", &
697 858 : pw_grid%dvol, " Volume (a.u.^3) :", pw_grid%vol
698 429 : IF (pw_grid%grid_span == HALFSPACE) THEN
699 289 : 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 28022 : n(1) = pw_grid%ngpts_cut_local
707 28022 : n(2) = pw_grid%ngpts_local
708 28022 : CALL pw_grid%para%group%sum(n(1:2))
709 84066 : n(3) = SUM(pw_grid%para%nyzray)
710 112088 : rv(:, 1) = REAL(n, KIND=dp)/REAL(pw_grid%para%group%num_pe, KIND=dp)
711 28022 : n(1) = pw_grid%ngpts_cut_local
712 28022 : n(2) = pw_grid%ngpts_local
713 28022 : CALL pw_grid%para%group%max(n(1:2))
714 84066 : n(3) = MAXVAL(pw_grid%para%nyzray)
715 112088 : rv(:, 2) = REAL(n, KIND=dp)
716 28022 : n(1) = pw_grid%ngpts_cut_local
717 28022 : n(2) = pw_grid%ngpts_local
718 28022 : CALL pw_grid%para%group%min(n(1:2))
719 84066 : n(3) = MINVAL(pw_grid%para%nyzray)
720 112088 : rv(:, 3) = REAL(n, KIND=dp)
721 :
722 28022 : IF (info > 0) THEN
723 : WRITE (info, '(/,A,T71,I10)') &
724 6289 : " PW_GRID| Information for grid number ", pw_grid%id_nr
725 6289 : IF (pw_grid%reference > 0) THEN
726 : WRITE (info, '(A,T71,I10)') &
727 4062 : " PW_GRID| Number of the reference grid ", pw_grid%reference
728 : END IF
729 : WRITE (info, '(A,T60,I10,A)') &
730 6289 : " PW_GRID| Grid distributed over ", pw_grid%para%group%num_pe, &
731 12578 : " processors"
732 : WRITE (info, '(A,T71,2I5)') &
733 6289 : " PW_GRID| Real space group dimensions ", pw_grid%para%group%num_pe_cart
734 6289 : IF (pw_grid%para%blocked) THEN
735 6 : WRITE (info, '(A,T78,A)') " PW_GRID| the grid is blocked: ", "YES"
736 : ELSE
737 6283 : WRITE (info, '(A,T78,A)') " PW_GRID| the grid is blocked: ", " NO"
738 : END IF
739 6289 : WRITE (info, '(" PW_GRID| Cutoff [a.u.]",T71,f10.1)') pw_grid%cutoff
740 6289 : IF (pw_grid%spherical) THEN
741 392 : WRITE (info, '(A,T78,A)') " PW_GRID| spherical cutoff: ", "YES"
742 392 : WRITE (info, '(A,T71,I10)') " PW_GRID| Grid points within cutoff", &
743 784 : pw_grid%ngpts_cut
744 : ELSE
745 5897 : WRITE (info, '(A,T78,A)') " PW_GRID| spherical cutoff: ", " NO"
746 : END IF
747 25156 : DO i = 1, 3
748 18867 : WRITE (info, '(A,I3,T30,2I8,T62,A,T71,I10)') " PW_GRID| Bounds ", &
749 18867 : i, pw_grid%bounds(1, I), pw_grid%bounds(2, I), &
750 44023 : "Points:", pw_grid%npts(I)
751 : END DO
752 : WRITE (info, '(A,G12.4,T50,A,T67,F14.4)') &
753 6289 : " PW_GRID| Volume element (a.u.^3)", &
754 12578 : pw_grid%dvol, " Volume (a.u.^3) :", pw_grid%vol
755 6289 : IF (pw_grid%grid_span == HALFSPACE) THEN
756 400 : WRITE (info, '(A,T72,A)') " PW_GRID| Grid span", "HALFSPACE"
757 : ELSE
758 5889 : WRITE (info, '(A,T72,A)') " PW_GRID| Grid span", "FULLSPACE"
759 : END IF
760 6289 : WRITE (info, '(A,T48,A)') " PW_GRID| Distribution", &
761 12578 : " Average Max Min"
762 6289 : WRITE (info, '(A,T45,F12.1,2I12)') " PW_GRID| G-Vectors", &
763 12578 : rv(1, 1), NINT(rv(1, 2)), NINT(rv(1, 3))
764 6289 : WRITE (info, '(A,T45,F12.1,2I12)') " PW_GRID| G-Rays ", &
765 12578 : rv(3, 1), NINT(rv(3, 2)), NINT(rv(3, 3))
766 6289 : WRITE (info, '(A,T45,F12.1,2I12)') " PW_GRID| Real Space Points", &
767 12578 : rv(2, 1), NINT(rv(2, 2)), NINT(rv(2, 3))
768 : END IF ! group head
769 : END IF ! local
770 :
771 33125 : 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 33147 : 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 33147 : INTEGER, ALLOCATABLE, DIMENSION(:) :: pemap
806 33147 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: yz_index
807 33147 : 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 33147 : CALL timeset(routineN, handle)
815 :
816 33147 : lby = pw_grid%bounds(1, 2)
817 33147 : lbz = pw_grid%bounds(1, 3)
818 :
819 132588 : pw_grid%ngpts = PRODUCT(INT(pw_grid%npts, KIND=int_8))
820 :
821 33147 : my_rs_dims = 0
822 33147 : IF (PRESENT(rs_dims)) THEN
823 32136 : my_rs_dims = rs_dims
824 : END IF
825 :
826 33147 : IF (PRESENT(blocked)) THEN
827 29766 : blocked_local = blocked
828 : ELSE
829 : blocked_local = do_pw_grid_blocked_free
830 : END IF
831 :
832 33147 : pw_grid%para%blocked = .FALSE.
833 :
834 33147 : IF (pw_grid%para%mode == PW_MODE_LOCAL) THEN
835 :
836 5103 : pw_grid%para%ray_distribution = .FALSE.
837 :
838 51030 : pw_grid%bounds_local = pw_grid%bounds
839 20412 : pw_grid%npts_local = pw_grid%npts
840 5103 : CPASSERT(pw_grid%ngpts_cut < HUGE(pw_grid%ngpts_cut_local))
841 5103 : pw_grid%ngpts_cut_local = INT(pw_grid%ngpts_cut)
842 5103 : CPASSERT(pw_grid%ngpts < HUGE(pw_grid%ngpts_local))
843 5103 : pw_grid%ngpts_local = INT(pw_grid%ngpts)
844 15309 : my_rs_dims = 1
845 5103 : CALL pw_grid%para%group%create(mp_comm, 2, my_rs_dims)
846 :
847 5103 : ALLOCATE (pw_grid%para%bo(2, 3, 0:0, 3))
848 20412 : DO i = 1, 3
849 61236 : pw_grid%para%bo(1, 1:3, 0, i) = 1
850 66339 : 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 28044 : nx = pw_grid%npts(1)
857 28044 : ny = pw_grid%npts(2)
858 28044 : nz = pw_grid%npts(3)
859 28044 : 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 39856 : IF (ANY(my_rs_dims <= 0)) THEN
871 66390 : IF (ALL(my_rs_dims <= 0)) THEN
872 22114 : 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 84132 : 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 28092 : IF (ALL(my_rs_dims == [1, np])) my_rs_dims = [0, 0]
885 :
886 : ! if [0,0] now, we can optimize it ourselves
887 72328 : IF (ALL(my_rs_dims == [0, 0])) THEN
888 : ! only small grids have a chance to be 2d distributed
889 22142 : 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 66426 : 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 26526 : IF (ALL(my_rs_dims == [np, 1])) THEN
917 : blocking = .FALSE.
918 : ELSE
919 0 : blocking = .TRUE.
920 : END IF
921 : CASE DEFAULT
922 28044 : CPABORT("")
923 : END SELECT
924 :
925 : !..create group for real space distribution
926 28044 : CALL pw_grid%para%group%create(mp_comm, 2, my_rs_dims)
927 :
928 28044 : 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 28022 : pw_grid%para%group%mepos_cart(1))
933 84066 : 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 28022 : pw_grid%para%group%mepos_cart(2))
936 84066 : pw_grid%bounds_local(:, 2) = lo + pw_grid%bounds(1, 2) - 1
937 84066 : pw_grid%bounds_local(:, 3) = pw_grid%bounds(:, 3)
938 : END IF
939 :
940 : pw_grid%npts_local(:) = pw_grid%bounds_local(2, :) &
941 112176 : - pw_grid%bounds_local(1, :) + 1
942 :
943 : !..the third distribution is needed for the second step in the FFT
944 112176 : ALLOCATE (pw_grid%para%bo(2, 3, 0:np - 1, 3))
945 84132 : rsd = pw_grid%para%group%num_pe_cart
946 :
947 28044 : 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 84066 : DO ip = 0, np - 1
975 56044 : CALL pw_grid%para%group%coords(ip, coor)
976 : ! distribution xyZ
977 168132 : pw_grid%para%bo(1:2, 1, ip, 1) = get_limit(nx, rsd(1), coor(1))
978 168132 : pw_grid%para%bo(1:2, 2, ip, 1) = get_limit(ny, rsd(2), coor(2))
979 56044 : pw_grid%para%bo(1, 3, ip, 1) = 1
980 56044 : pw_grid%para%bo(2, 3, ip, 1) = nz
981 : ! distribution xYz
982 168132 : pw_grid%para%bo(1:2, 1, ip, 2) = get_limit(nx, rsd(1), coor(1))
983 56044 : pw_grid%para%bo(1, 2, ip, 2) = 1
984 56044 : pw_grid%para%bo(2, 2, ip, 2) = ny
985 168132 : pw_grid%para%bo(1:2, 3, ip, 2) = get_limit(nz, rsd(2), coor(2))
986 : ! distribution Xyz
987 56044 : pw_grid%para%bo(1, 1, ip, 3) = 1
988 56044 : pw_grid%para%bo(2, 1, ip, 3) = nx
989 168132 : pw_grid%para%bo(1:2, 2, ip, 3) = get_limit(ny, rsd(1), coor(1))
990 196154 : 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 28044 : pw_grid%ngpts_cut_local = 0
995 :
996 84132 : ALLOCATE (pw_grid%para%nyzray(0:np - 1))
997 :
998 112176 : ALLOCATE (pw_grid%para%yzq(ny, nz))
999 :
1000 : IF (pw_grid%spherical .OR. pw_grid%grid_span == HALFSPACE &
1001 28044 : .OR. .NOT. blocking) THEN
1002 :
1003 28032 : pw_grid%para%ray_distribution = .TRUE.
1004 :
1005 31147498 : pw_grid%para%yzq = -1
1006 28032 : IF (PRESENT(ref_grid)) THEN
1007 : ! tag all vectors from the reference grid
1008 17948 : 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 28032 : i1 = SIZE(yz_mask, 1)
1016 28032 : i2 = SIZE(yz_mask, 2)
1017 84096 : ALLOCATE (yz_index(2, i1*i2))
1018 28032 : CALL order_mask(yz_mask, yz_index)
1019 30398592 : DO i = 1, i1*i2
1020 30370560 : lo(1) = yz_index(1, i)
1021 30370560 : lo(2) = yz_index(2, i)
1022 30370560 : IF (lo(1)*lo(2) == 0) CYCLE
1023 18682542 : gmax = yz_mask(lo(1), lo(2))
1024 18682542 : IF (gmax == 0) CYCLE
1025 18682542 : yz_mask(lo(1), lo(2)) = 0
1026 18682542 : ip = MOD(i - 1, 2*np)
1027 18682542 : IF (ip > np - 1) ip = 2*np - ip - 1
1028 18682542 : IF (ip == pw_grid%para%group%mepos) THEN
1029 9341271 : pw_grid%ngpts_cut_local = pw_grid%ngpts_cut_local + gmax
1030 : END IF
1031 18682542 : pw_grid%para%yzq(lo(1), lo(2)) = ip
1032 18710574 : IF (pw_grid%grid_span == HALFSPACE) THEN
1033 1227096 : m = -lo(1) - 2*lby + 2
1034 1227096 : n = -lo(2) - 2*lbz + 2
1035 1227096 : pw_grid%para%yzq(m, n) = ip
1036 1227096 : yz_mask(m, n) = 0
1037 : END IF
1038 : END DO
1039 :
1040 28032 : DEALLOCATE (yz_index)
1041 :
1042 : ! Count the total number of rays on each processor
1043 84096 : pw_grid%para%nyzray = 0
1044 776938 : DO i = 1, nz
1045 31147498 : DO j = 1, ny
1046 30370560 : ip = pw_grid%para%yzq(j, i)
1047 30370560 : IF (ip >= 0) pw_grid%para%nyzray(ip) = &
1048 29917018 : pw_grid%para%nyzray(ip) + 1
1049 : END DO
1050 : END DO
1051 :
1052 : ! Allocate mapping array (y:z, nray, nproc)
1053 84096 : ns = MAXVAL(pw_grid%para%nyzray(0:np - 1))
1054 112128 : ALLOCATE (pw_grid%para%yzp(2, ns, 0:np - 1))
1055 :
1056 : ! Fill mapping array, recalculate nyzray for convenience
1057 84096 : pw_grid%para%nyzray = 0
1058 776938 : DO i = 1, nz
1059 31147498 : DO j = 1, ny
1060 30370560 : ip = pw_grid%para%yzq(j, i)
1061 31119466 : IF (ip >= 0) THEN
1062 : pw_grid%para%nyzray(ip) = &
1063 29168112 : pw_grid%para%nyzray(ip) + 1
1064 29168112 : ns = pw_grid%para%nyzray(ip)
1065 29168112 : pw_grid%para%yzp(1, ns, ip) = j
1066 29168112 : pw_grid%para%yzp(2, ns, ip) = i
1067 29168112 : IF (ip == pw_grid%para%group%mepos) THEN
1068 14584056 : pw_grid%para%yzq(j, i) = ns
1069 : ELSE
1070 14584056 : pw_grid%para%yzq(j, i) = -1
1071 : END IF
1072 : ELSE
1073 1202448 : pw_grid%para%yzq(j, i) = -2
1074 : END IF
1075 : END DO
1076 : END DO
1077 :
1078 112128 : 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 33147 : IF (pw_grid%para%mode .EQ. PW_MODE_DISTRIBUTED) THEN
1155 84132 : ALLOCATE (pw_grid%para%pos_of_x(pw_grid%bounds(1, 1):pw_grid%bounds(2, 1)))
1156 768604 : pw_grid%para%pos_of_x = 0
1157 398324 : pw_grid%para%pos_of_x(pw_grid%bounds_local(1, 1):pw_grid%bounds_local(2, 1)) = pw_grid%para%group%mepos
1158 28044 : CALL pw_grid%para%group%sum(pw_grid%para%pos_of_x)
1159 : ELSE
1160 : ! this should not be needed
1161 15309 : ALLOCATE (pw_grid%para%pos_of_x(pw_grid%bounds(1, 1):pw_grid%bounds(2, 1)))
1162 86334 : pw_grid%para%pos_of_x = 0
1163 : END IF
1164 :
1165 33147 : CALL timestop(handle)
1166 :
1167 33147 : 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 17948 : 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 17948 : ny = ref_grid%npts(2)
1187 17948 : nz = ref_grid%npts(3)
1188 17948 : lby = pw_grid%bounds(1, 2)
1189 17948 : lbz = pw_grid%bounds(1, 3)
1190 17948 : uby = pw_grid%bounds(2, 2)
1191 17948 : ubz = pw_grid%bounds(2, 3)
1192 17948 : my = SIZE(yz_mask, 1)
1193 17948 : mz = SIZE(yz_mask, 2)
1194 :
1195 : ! loop over all processors and all g vectors yz lines on this processor
1196 53844 : DO ip = 0, ref_grid%para%group%num_pe - 1
1197 47476162 : 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 47422318 : y = ref_grid%para%yzp(1, ig, ip) - 1
1201 47422318 : IF (y >= ny/2) y = y - ny
1202 47422318 : z = ref_grid%para%yzp(2, ig, ip) - 1
1203 47422318 : IF (z >= nz/2) z = z - nz
1204 : ! check if this is inside the realm of the new grid
1205 47422318 : IF (y < lby .OR. y > uby .OR. z < lbz .OR. z > ubz) CYCLE
1206 : ! go to shifted coordinates
1207 9260734 : y = y - lby + 1
1208 9260734 : z = z - lbz + 1
1209 : ! this tag is outside the cutoff range of the new grid
1210 9260734 : 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 9260734 : gmax = yz_mask(y, z)
1224 9260734 : IF (gmax == 0) CYCLE
1225 9260734 : yz_mask(y, z) = 0
1226 9260734 : pw_grid%para%yzq(y, z) = ip
1227 : END IF
1228 9296630 : IF (ip == pw_grid%para%group%mepos) THEN
1229 4630367 : pw_grid%ngpts_cut_local = pw_grid%ngpts_cut_local + gmax
1230 : END IF
1231 : END DO
1232 : END DO
1233 :
1234 17948 : END SUBROUTINE pre_tag
1235 :
1236 : ! **************************************************************************************************
1237 : !> \brief ...
1238 : !> \param yz_mask ...
1239 : !> \param yz_index ...
1240 : ! **************************************************************************************************
1241 28032 : 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 28032 : i1 = SIZE(yz_mask, 1)
1255 28032 : i2 = SIZE(yz_mask, 2)
1256 91139712 : yz_index = 0
1257 :
1258 28032 : icount = 1
1259 28032 : ic = i1/2
1260 28032 : jc = i2/2
1261 28032 : ii = ic
1262 28032 : jj = jc
1263 28032 : IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1264 28032 : IF (yz_mask(ii, jj) /= 0) THEN
1265 10084 : yz_index(1, icount) = ii
1266 10084 : yz_index(2, icount) = jj
1267 10084 : icount = icount + 1
1268 : END IF
1269 : END IF
1270 428462 : DO im = 1, MAX(ic + 1, jc + 1)
1271 400430 : ii = ic - im
1272 10362956 : DO jj = jc - im, jc + im
1273 10362956 : IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1274 7297796 : IF (yz_mask(ii, jj) /= 0) THEN
1275 4549730 : yz_index(1, icount) = ii
1276 4549730 : yz_index(2, icount) = jj
1277 4549730 : icount = icount + 1
1278 : END IF
1279 : END IF
1280 : END DO
1281 400430 : ii = ic + im
1282 10362956 : DO jj = jc - im, jc + im
1283 10362956 : IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1284 8231796 : IF (yz_mask(ii, jj) /= 0) THEN
1285 4946532 : yz_index(1, icount) = ii
1286 4946532 : yz_index(2, icount) = jj
1287 4946532 : icount = icount + 1
1288 : END IF
1289 : END IF
1290 : END DO
1291 400430 : jj = jc - im
1292 9562096 : DO ii = ic - im + 1, ic + im - 1
1293 9562096 : IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1294 6958780 : IF (yz_mask(ii, jj) /= 0) THEN
1295 4675098 : yz_index(1, icount) = ii
1296 4675098 : yz_index(2, icount) = jj
1297 4675098 : icount = icount + 1
1298 : END IF
1299 : END IF
1300 : END DO
1301 9562096 : jj = jc + im
1302 9590128 : DO ii = ic - im + 1, ic + im - 1
1303 9562096 : IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1304 7854156 : IF (yz_mask(ii, jj) /= 0) THEN
1305 4501098 : yz_index(1, icount) = ii
1306 4501098 : yz_index(2, icount) = jj
1307 4501098 : icount = icount + 1
1308 : END IF
1309 : END IF
1310 : END DO
1311 : END DO
1312 :
1313 28032 : 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 1688121342 : 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 1688121342 : + 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 1688121342 : + 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 1688121342 : + REAL(n, dp)*h_inv(3, 3)
1343 :
1344 : ! enforce strict zero-ness in this case (compiler optimization)
1345 1688121342 : IF (l == 0 .AND. m == 0 .AND. n == 0) THEN
1346 38250 : length_x = 0
1347 38250 : length_y = 0
1348 38250 : length_z = 0
1349 : END IF
1350 :
1351 1688121342 : length_x = length_x*twopi
1352 1688121342 : length_y = length_y*twopi
1353 1688121342 : length_z = length_z*twopi
1354 :
1355 1688121342 : length = length_x**2 + length_y**2 + length_z**2
1356 :
1357 1688121342 : 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 33147 : 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 33147 : IF (pw_grid%grid_span == HALFSPACE) THEN
1387 : n_upperlimit = 0
1388 29736 : ELSE IF (pw_grid%grid_span == FULLSPACE) THEN
1389 29736 : 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 33147 : gpt = 0
1396 33147 : IF (pw_grid%para%mode == PW_MODE_LOCAL) THEN
1397 5103 : nlim(1) = bounds(1, 3)
1398 5103 : nlim(2) = n_upperlimit
1399 28044 : ELSE IF (pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
1400 28044 : n = n_upperlimit - bounds(1, 3) + 1
1401 28044 : nlim = get_limit(n, mp_comm%num_pe, mp_comm%mepos)
1402 84132 : nlim = nlim + bounds(1, 3) - 1
1403 : ELSE
1404 0 : CPABORT("para % mode not specified")
1405 : END IF
1406 :
1407 33435067 : yz_mask = 0
1408 486051 : DO n = nlim(1), nlim(2)
1409 419757 : nn = n - bounds(1, 3) + 1
1410 16267353 : DO m = bounds(1, 2), bounds(2, 2)
1411 15814449 : mm = m - bounds(1, 2) + 1
1412 869476559 : DO l = bounds(1, 1), bounds(2, 1)
1413 853242353 : IF (pw_grid%grid_span == HALFSPACE .AND. n == 0) THEN
1414 2965529 : IF ((m == 0 .AND. l > 0) .OR. (m > 0)) CYCLE
1415 : END IF
1416 :
1417 851649296 : CALL pw_vec_length(h_inv, length_x, length_y, length_z, length, l, m, n)
1418 :
1419 867463745 : IF (0.5_dp*length <= cutoff) THEN
1420 811001030 : gpt = gpt + 1
1421 811001030 : 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 33147 : IF (pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
1431 28044 : CALL mp_comm%sum(gpt)
1432 62276584 : CALL mp_comm%sum(yz_mask)
1433 : END IF
1434 33147 : pw_grid%ngpts_cut = gpt
1435 :
1436 33147 : 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 33147 : 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 33147 : CALL timeset(routineN, handle)
1463 :
1464 331470 : bounds = pw_grid%bounds
1465 33147 : lby = pw_grid%bounds(1, 2)
1466 33147 : lbz = pw_grid%bounds(1, 3)
1467 :
1468 33147 : IF (pw_grid%grid_span == HALFSPACE) THEN
1469 : n_upperlimit = 0
1470 29736 : ELSE IF (pw_grid%grid_span == FULLSPACE) THEN
1471 29736 : 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 33147 : IF (pw_grid%para%mode == PW_MODE_LOCAL) THEN
1478 5103 : gpt = 0
1479 70085 : DO n = bounds(1, 3), n_upperlimit
1480 1579939 : DO m = bounds(1, 2), bounds(2, 2)
1481 59220684 : DO l = bounds(1, 1), bounds(2, 1)
1482 57645848 : IF (pw_grid%grid_span == HALFSPACE .AND. n == 0) THEN
1483 1157021 : IF ((m == 0 .AND. l > 0) .OR. (m > 0)) CYCLE
1484 : END IF
1485 :
1486 56946134 : CALL pw_vec_length(h_inv, length_x, length_y, length_z, length, l, m, n)
1487 :
1488 58455988 : IF (0.5_dp*length <= cutoff) THEN
1489 43347699 : gpt = gpt + 1
1490 43347699 : pw_grid%g(1, gpt) = length_x
1491 43347699 : pw_grid%g(2, gpt) = length_y
1492 43347699 : pw_grid%g(3, gpt) = length_z
1493 43347699 : pw_grid%gsq(gpt) = length
1494 43347699 : pw_grid%g_hat(1, gpt) = l
1495 43347699 : pw_grid%g_hat(2, gpt) = m
1496 43347699 : 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 28044 : IF (pw_grid%para%ray_distribution) THEN
1506 :
1507 28032 : gpt = 0
1508 28032 : ip = pw_grid%para%group%mepos
1509 14612088 : DO i = 1, pw_grid%para%nyzray(ip)
1510 14584056 : n = pw_grid%para%yzp(2, i, ip) + lbz - 1
1511 14584056 : m = pw_grid%para%yzp(1, i, ip) + lby - 1
1512 14584056 : IF (n > n_upperlimit) CYCLE
1513 794314035 : DO l = bounds(1, 1), bounds(2, 1)
1514 780296779 : IF (pw_grid%grid_span == HALFSPACE .AND. n == 0) THEN
1515 1629336 : IF ((m == 0 .AND. l > 0) .OR. (m > 0)) CYCLE
1516 : END IF
1517 :
1518 779482991 : CALL pw_vec_length(h_inv, length_x, length_y, length_z, length, l, m, n)
1519 :
1520 794067047 : IF (0.5_dp*length <= cutoff) THEN
1521 767610410 : gpt = gpt + 1
1522 767610410 : pw_grid%g(1, gpt) = length_x
1523 767610410 : pw_grid%g(2, gpt) = length_y
1524 767610410 : pw_grid%g(3, gpt) = length_z
1525 767610410 : pw_grid%gsq(gpt) = length
1526 767610410 : pw_grid%g_hat(1, gpt) = l
1527 767610410 : pw_grid%g_hat(2, gpt) = m
1528 767610410 : 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 33147 : CPASSERT(pw_grid%ngpts_cut_local == gpt)
1581 33147 : IF (pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
1582 28044 : gpt_global = gpt
1583 28044 : CALL pw_grid%para%group%sum(gpt_global)
1584 28044 : CPASSERT(pw_grid%ngpts_cut == gpt_global)
1585 : END IF
1586 :
1587 33147 : pw_grid%have_g0 = .FALSE.
1588 33147 : pw_grid%first_gne0 = 1
1589 611524609 : DO gpt = 1, pw_grid%ngpts_cut_local
1590 623427317 : IF (ALL(pw_grid%g_hat(:, gpt) == 0)) THEN
1591 19125 : pw_grid%have_g0 = .TRUE.
1592 19125 : pw_grid%first_gne0 = 2
1593 19125 : EXIT
1594 : END IF
1595 : END DO
1596 :
1597 : CALL pw_grid_set_maps(pw_grid%grid_span, pw_grid%g_hat, &
1598 33147 : pw_grid%mapl, pw_grid%mapm, pw_grid%mapn, pw_grid%npts)
1599 :
1600 33147 : CALL timestop(handle)
1601 :
1602 33147 : 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 33147 : 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 33147 : ng = SIZE(g_hat, 2)
1629 :
1630 811034177 : DO gpt = 1, ng
1631 811001030 : l = g_hat(1, gpt)
1632 811001030 : m = g_hat(2, gpt)
1633 811001030 : n = g_hat(3, gpt)
1634 811001030 : IF (l < 0) THEN
1635 403206492 : mapl%pos(l) = l + npts(1)
1636 : ELSE
1637 407794538 : mapl%pos(l) = l
1638 : END IF
1639 811001030 : IF (m < 0) THEN
1640 403628183 : mapm%pos(m) = m + npts(2)
1641 : ELSE
1642 407372847 : mapm%pos(m) = m
1643 : END IF
1644 811001030 : IF (n < 0) THEN
1645 419365092 : mapn%pos(n) = n + npts(3)
1646 : ELSE
1647 391635938 : mapn%pos(n) = n
1648 : END IF
1649 :
1650 : ! Generating the maps to the full 3-d space
1651 :
1652 811034177 : IF (grid_span == HALFSPACE) THEN
1653 :
1654 33226993 : IF (l <= 0) THEN
1655 17102585 : mapl%neg(l) = -l
1656 : ELSE
1657 16124408 : mapl%neg(l) = npts(1) - l
1658 : END IF
1659 33226993 : IF (m <= 0) THEN
1660 17538924 : mapm%neg(m) = -m
1661 : ELSE
1662 15688069 : mapm%neg(m) = npts(2) - m
1663 : END IF
1664 33226993 : IF (n <= 0) THEN
1665 33226993 : 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 33147 : 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 33147 : 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 33147 : CALL timeset(routineN, handle)
1707 :
1708 99441 : ALLOCATE (pw_grid%g(3, ng))
1709 99441 : ALLOCATE (pw_grid%gsq(ng))
1710 99441 : ALLOCATE (pw_grid%g_hat(3, ng))
1711 :
1712 : nmaps = 1
1713 33147 : 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 33147 : ALLOCATE (pw_grid%g_hatmap(1, 1))
1723 : #endif
1724 :
1725 33147 : IF (pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
1726 : ALLOCATE (pw_grid%grays(pw_grid%npts(1), &
1727 112176 : pw_grid%para%nyzray(pw_grid%para%group%mepos)))
1728 : END IF
1729 :
1730 99441 : ALLOCATE (pw_grid%mapl%pos(bounds(1, 1):bounds(2, 1)))
1731 99441 : ALLOCATE (pw_grid%mapl%neg(bounds(1, 1):bounds(2, 1)))
1732 99441 : ALLOCATE (pw_grid%mapm%pos(bounds(1, 2):bounds(2, 2)))
1733 99441 : ALLOCATE (pw_grid%mapm%neg(bounds(1, 2):bounds(2, 2)))
1734 99441 : ALLOCATE (pw_grid%mapn%pos(bounds(1, 3):bounds(2, 3)))
1735 99441 : ALLOCATE (pw_grid%mapn%neg(bounds(1, 3):bounds(2, 3)))
1736 :
1737 33147 : CALL timestop(handle)
1738 :
1739 33147 : 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 33147 : 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 33147 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: int_tmp
1771 : LOGICAL :: g_found
1772 : REAL(KIND=dp) :: gig, gigr
1773 33147 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: real_tmp
1774 :
1775 33147 : CALL timeset(routineN, handle)
1776 :
1777 33147 : ng = SIZE(pw_grid%gsq)
1778 99441 : ALLOCATE (idx(ng))
1779 :
1780 : ! grids are (locally) ordered by length of G-vectors
1781 33147 : CALL sort(pw_grid%gsq, ng, idx)
1782 : ! within shells order wrt x,y,z
1783 33147 : CALL sort_shells(pw_grid%gsq, pw_grid%g_hat, idx)
1784 :
1785 99441 : ALLOCATE (real_tmp(3, ng))
1786 811034177 : DO i = 1, ng
1787 811001030 : real_tmp(1, i) = pw_grid%g(1, idx(i))
1788 811001030 : real_tmp(2, i) = pw_grid%g(2, idx(i))
1789 811034177 : real_tmp(3, i) = pw_grid%g(3, idx(i))
1790 : END DO
1791 811034177 : DO i = 1, ng
1792 811001030 : pw_grid%g(1, i) = real_tmp(1, i)
1793 811001030 : pw_grid%g(2, i) = real_tmp(2, i)
1794 811034177 : pw_grid%g(3, i) = real_tmp(3, i)
1795 : END DO
1796 33147 : DEALLOCATE (real_tmp)
1797 :
1798 99441 : ALLOCATE (int_tmp(3, ng))
1799 811034177 : DO i = 1, ng
1800 811001030 : int_tmp(1, i) = pw_grid%g_hat(1, idx(i))
1801 811001030 : int_tmp(2, i) = pw_grid%g_hat(2, idx(i))
1802 811034177 : int_tmp(3, i) = pw_grid%g_hat(3, idx(i))
1803 : END DO
1804 811034177 : DO i = 1, ng
1805 811001030 : pw_grid%g_hat(1, i) = int_tmp(1, i)
1806 811001030 : pw_grid%g_hat(2, i) = int_tmp(2, i)
1807 811034177 : pw_grid%g_hat(3, i) = int_tmp(3, i)
1808 : END DO
1809 33147 : DEALLOCATE (int_tmp)
1810 :
1811 33147 : DEALLOCATE (idx)
1812 :
1813 : ! check if ordering is compatible to reference grid
1814 33147 : IF (PRESENT(ref_grid)) THEN
1815 20912 : ngr = SIZE(ref_grid%gsq)
1816 20912 : ngr = MIN(ng, ngr)
1817 20912 : 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 62736 : ALLOCATE (pw_grid%gidx(1:ngr))
1824 176899881 : pw_grid%gidx = 0
1825 : ! first try as many trivial associations as possible
1826 20912 : it = 0
1827 94002427 : DO ig = 1, ngr
1828 375949780 : IF (.NOT. ALL(pw_grid%g_hat(1:3, ig) &
1829 : == ref_grid%g_hat(1:3, ig))) EXIT
1830 93981515 : pw_grid%gidx(ig) = ig
1831 94002427 : it = ig
1832 : END DO
1833 : ! now for the ones that could not be done
1834 20912 : IF (ng == ngr) THEN
1835 : ! for the case pw_grid <= ref_grid
1836 20912 : is = it
1837 82918366 : DO ig = it + 1, ngr
1838 82897454 : gig = pw_grid%gsq(ig)
1839 82897454 : gigr = MAX(1._dp, gig)
1840 82897454 : g_found = .FALSE.
1841 416340475 : DO ih = is + 1, SIZE(ref_grid%gsq)
1842 416340475 : IF (ABS(gig - ref_grid%gsq(ih))/gigr > 1.e-12_dp) CYCLE
1843 : g_found = .TRUE.
1844 333443021 : EXIT
1845 : END DO
1846 82897454 : 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 82897454 : ip = ih - 1
1851 5698733083 : DO ih = ip + 1, SIZE(ref_grid%gsq)
1852 5698733083 : IF (ABS(gig - ref_grid%gsq(ih))/gigr > 1.e-12_dp) CYCLE
1853 5698733083 : IF (pw_grid%g_hat(1, ig) /= ref_grid%g_hat(1, ih)) CYCLE
1854 316146545 : IF (pw_grid%g_hat(2, ig) /= ref_grid%g_hat(2, ih)) CYCLE
1855 102588510 : IF (pw_grid%g_hat(3, ig) /= ref_grid%g_hat(3, ih)) CYCLE
1856 82897454 : pw_grid%gidx(ig) = ih
1857 5698733083 : EXIT
1858 : END DO
1859 82897454 : 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 82918366 : 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 176899881 : 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 33147 : IF (pw_grid%have_g0) THEN
1924 19125 : 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 33147 : CALL timestop(handle)
1930 :
1931 33147 : END SUBROUTINE pw_grid_sort
1932 :
1933 : ! **************************************************************************************************
1934 : !> \brief ...
1935 : !> \param gsq ...
1936 : !> \param g_hat ...
1937 : !> \param idx ...
1938 : ! **************************************************************************************************
1939 33147 : 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 33147 : 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 33147 : ng = SIZE(gsq)
1958 33147 : s_begin = -1.0_dp
1959 33147 : s1 = 0
1960 33147 : s2 = 0
1961 : ig = 1
1962 811034177 : DO ig = 1, ng
1963 811034177 : IF (ABS(gsq(ig) - s_begin) < small) THEN
1964 770549400 : s2 = ig
1965 : ELSE
1966 40451630 : CALL redist(g_hat, idx, s1, s2)
1967 40451630 : s_begin = gsq(ig)
1968 40451630 : s1 = ig
1969 40451630 : s2 = ig
1970 : END IF
1971 : END DO
1972 33147 : CALL redist(g_hat, idx, s1, s2)
1973 :
1974 33147 : CALL timestop(handle)
1975 :
1976 33147 : END SUBROUTINE sort_shells
1977 :
1978 : ! **************************************************************************************************
1979 : !> \brief ...
1980 : !> \param g_hat ...
1981 : !> \param idx ...
1982 : !> \param s1 ...
1983 : !> \param s2 ...
1984 : ! **************************************************************************************************
1985 40484777 : 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 40484777 : INTEGER, ALLOCATABLE, DIMENSION(:) :: indl
1994 40484777 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: slen
1995 :
1996 40484777 : IF (s2 <= s1) RETURN
1997 38704670 : ns = s2 - s1 + 1
1998 116114010 : ALLOCATE (indl(ns))
1999 116114010 : ALLOCATE (slen(ns))
2000 :
2001 847958740 : DO i = s1, s2
2002 809254070 : ii = idx(i)
2003 809254070 : n1 = g_hat(1, ii)
2004 809254070 : n2 = g_hat(2, ii)
2005 809254070 : n3 = g_hat(3, ii)
2006 : slen(i - s1 + 1) = 1000.0_dp*REAL(n1, dp) + &
2007 847958740 : REAL(n2, dp) + 0.001_dp*REAL(n3, dp)
2008 : END DO
2009 38704670 : CALL sort(slen, ns, indl)
2010 847958740 : DO i = 1, ns
2011 809254070 : ii = indl(i) + s1 - 1
2012 847958740 : indl(i) = idx(ii)
2013 : END DO
2014 847958740 : idx(s1:s2) = indl(1:ns)
2015 :
2016 38704670 : DEALLOCATE (indl)
2017 38704670 : 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 33147 : 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 33147 : IF (pw_grid%para%mode == PW_MODE_LOCAL) RETURN
2040 :
2041 28044 : 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 31152314 : yz = 0
2047 87696588 : pw_grid%para%yzp = 0
2048 31152314 : pw_grid%para%yzq = 0
2049 :
2050 767681375 : DO gpt = 1, SIZE(pw_grid%gsq)
2051 767653331 : m = posm(pw_grid%g_hat(2, gpt)) + 1
2052 767653331 : n = posn(pw_grid%g_hat(3, gpt)) + 1
2053 767681375 : yz(m, n) = yz(m, n) + 1
2054 : END DO
2055 28044 : IF (pw_grid%grid_span == HALFSPACE) THEN
2056 20665962 : DO gpt = 1, SIZE(pw_grid%gsq)
2057 20663702 : m = negm(pw_grid%g_hat(2, gpt)) + 1
2058 20663702 : n = negn(pw_grid%g_hat(3, gpt)) + 1
2059 20665962 : yz(m, n) = yz(m, n) + 1
2060 : END DO
2061 : END IF
2062 :
2063 28044 : ip = pw_grid%para%group%mepos
2064 28044 : is = 0
2065 805236 : DO i = 1, nz
2066 31152314 : DO j = 1, ny
2067 31124270 : IF (yz(j, i) > 0) THEN
2068 14586337 : is = is + 1
2069 14586337 : pw_grid%para%yzp(1, is, ip) = j
2070 14586337 : pw_grid%para%yzp(2, is, ip) = i
2071 14586337 : pw_grid%para%yzq(j, i) = is
2072 : END IF
2073 : END DO
2074 : END DO
2075 : END ASSOCIATE
2076 :
2077 28044 : CPASSERT(is == pw_grid%para%nyzray(ip))
2078 28044 : CALL pw_grid%para%group%sum(pw_grid%para%yzp)
2079 :
2080 28044 : CALL timestop(handle)
2081 :
2082 28044 : 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 11114 : 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 11114 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: g
2106 :
2107 11114 : cell_deth = ABS(det_3x3(cell_hmat))
2108 11114 : 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 11114 : cell_h_inv = inv_3x3(cell_hmat)
2114 :
2115 11114 : g => pw_grid%g
2116 :
2117 11114 : CALL cell2grid(cell_hmat, cell_h_inv, cell_deth, pw_grid)
2118 :
2119 105197171 : DO gpt = 1, SIZE(g, 2)
2120 :
2121 105186057 : l = twopi*REAL(pw_grid%g_hat(1, gpt), KIND=dp)
2122 105186057 : m = twopi*REAL(pw_grid%g_hat(2, gpt), KIND=dp)
2123 105186057 : n = twopi*REAL(pw_grid%g_hat(3, gpt), KIND=dp)
2124 :
2125 105186057 : g(1, gpt) = l*cell_h_inv(1, 1) + m*cell_h_inv(2, 1) + n*cell_h_inv(3, 1)
2126 105186057 : g(2, gpt) = l*cell_h_inv(1, 2) + m*cell_h_inv(2, 2) + n*cell_h_inv(3, 2)
2127 105186057 : 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 105197171 : + g(3, gpt)*g(3, gpt)
2132 :
2133 : END DO
2134 :
2135 11114 : 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 7222410 : SUBROUTINE pw_grid_retain(pw_grid)
2147 : TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
2148 :
2149 7222410 : CPASSERT(pw_grid%ref_count > 0)
2150 7222410 : pw_grid%ref_count = pw_grid%ref_count + 1
2151 7222410 : 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 14864920 : 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 14864920 : IF (ASSOCIATED(pw_grid)) THEN
2172 7315631 : CPASSERT(pw_grid%ref_count > 0)
2173 7315631 : pw_grid%ref_count = pw_grid%ref_count - 1
2174 7315631 : IF (pw_grid%ref_count == 0) THEN
2175 93221 : IF (ASSOCIATED(pw_grid%gidx)) THEN
2176 20912 : DEALLOCATE (pw_grid%gidx)
2177 : END IF
2178 93221 : IF (ASSOCIATED(pw_grid%g)) THEN
2179 33147 : DEALLOCATE (pw_grid%g)
2180 : END IF
2181 93221 : IF (ASSOCIATED(pw_grid%gsq)) THEN
2182 33147 : DEALLOCATE (pw_grid%gsq)
2183 : END IF
2184 93221 : IF (ALLOCATED(pw_grid%g_hat)) THEN
2185 33147 : DEALLOCATE (pw_grid%g_hat)
2186 : END IF
2187 93221 : 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 33147 : DEALLOCATE (pw_grid%g_hatmap)
2194 : #endif
2195 : END IF
2196 93221 : IF (ASSOCIATED(pw_grid%grays)) THEN
2197 28044 : DEALLOCATE (pw_grid%grays)
2198 : END IF
2199 93221 : IF (ALLOCATED(pw_grid%mapl%pos)) THEN
2200 33147 : DEALLOCATE (pw_grid%mapl%pos)
2201 : END IF
2202 93221 : IF (ALLOCATED(pw_grid%mapm%pos)) THEN
2203 33147 : DEALLOCATE (pw_grid%mapm%pos)
2204 : END IF
2205 93221 : IF (ALLOCATED(pw_grid%mapn%pos)) THEN
2206 33147 : DEALLOCATE (pw_grid%mapn%pos)
2207 : END IF
2208 93221 : IF (ALLOCATED(pw_grid%mapl%neg)) THEN
2209 33147 : DEALLOCATE (pw_grid%mapl%neg)
2210 : END IF
2211 93221 : IF (ALLOCATED(pw_grid%mapm%neg)) THEN
2212 33147 : DEALLOCATE (pw_grid%mapm%neg)
2213 : END IF
2214 93221 : IF (ALLOCATED(pw_grid%mapn%neg)) THEN
2215 33147 : DEALLOCATE (pw_grid%mapn%neg)
2216 : END IF
2217 93221 : IF (ALLOCATED(pw_grid%para%bo)) THEN
2218 33147 : DEALLOCATE (pw_grid%para%bo)
2219 : END IF
2220 93221 : IF (pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
2221 28052 : IF (ALLOCATED(pw_grid%para%yzp)) THEN
2222 28044 : DEALLOCATE (pw_grid%para%yzp)
2223 : END IF
2224 28052 : IF (ALLOCATED(pw_grid%para%yzq)) THEN
2225 28044 : DEALLOCATE (pw_grid%para%yzq)
2226 : END IF
2227 28052 : IF (ALLOCATED(pw_grid%para%nyzray)) THEN
2228 28044 : DEALLOCATE (pw_grid%para%nyzray)
2229 : END IF
2230 : END IF
2231 : ! also release groups
2232 93221 : CALL pw_grid%para%group%free()
2233 93221 : IF (ALLOCATED(pw_grid%para%pos_of_x)) THEN
2234 34155 : DEALLOCATE (pw_grid%para%pos_of_x)
2235 : END IF
2236 :
2237 93221 : IF (ASSOCIATED(pw_grid)) THEN
2238 93221 : DEALLOCATE (pw_grid)
2239 : END IF
2240 : END IF
2241 : END IF
2242 14864920 : NULLIFY (pw_grid)
2243 14864920 : END SUBROUTINE pw_grid_release
2244 :
2245 : END MODULE pw_grids
|
