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