(git:374b731)
Loading...
Searching...
No Matches
pw_grids.F
Go to the documentation of this file.
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! **************************************************************************************************
36MODULE 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
81CONTAINS
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 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 my_local = .false.
102 IF (PRESENT(local)) my_local = local
103 cpassert(.NOT. ASSOCIATED(pw_grid))
104 ALLOCATE (pw_grid)
105 pw_grid%bounds = 0
106 pw_grid%cutoff = 0.0_dp
107 pw_grid%grid_span = fullspace
108 pw_grid%para%mode = pw_mode_local
109 pw_grid%para%rs_dims = 0
110 pw_grid%reference = 0
111 pw_grid%ref_count = 1
112 NULLIFY (pw_grid%g)
113 NULLIFY (pw_grid%gsq)
114 NULLIFY (pw_grid%g_hatmap)
115 NULLIFY (pw_grid%gidx)
116 NULLIFY (pw_grid%grays)
117
118 ! assign a unique tag to this grid
119 grid_tag = grid_tag + 1
120 pw_grid%id_nr = grid_tag
121
122 ! parallel info
123 CALL pw_grid%para%group%from_dup(pe_group)
124 pw_grid%para%group_size = pw_grid%para%group%num_pe
125 pw_grid%para%my_pos = pw_grid%para%group%mepos
126 pw_grid%para%group_head_id = 0
127 pw_grid%para%group_head = &
128 (pw_grid%para%group_head_id == pw_grid%para%my_pos)
129 IF (pw_grid%para%group_size > 1 .AND. (.NOT. my_local)) THEN
130 pw_grid%para%mode = pw_mode_distributed
131 ELSE
132 pw_grid%para%mode = pw_mode_local
133 END IF
134
135 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 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 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 equal = .false.
159 END IF
160
161 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 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 cpassert(pw_grid%ref_count > 0)
197
198 IF (PRESENT(id_nr)) id_nr = pw_grid%id_nr
199 IF (PRESENT(mode)) mode = pw_grid%para%mode
200 IF (PRESENT(vol)) vol = pw_grid%vol
201 IF (PRESENT(dvol)) dvol = pw_grid%dvol
202 IF (PRESENT(npts)) npts(1:3) = pw_grid%npts(1:3)
203 IF (PRESENT(ngpts)) ngpts = pw_grid%ngpts
204 IF (PRESENT(ngpts_cut)) ngpts_cut = pw_grid%ngpts_cut
205 IF (PRESENT(dr)) dr = pw_grid%dr
206 IF (PRESENT(cutoff)) cutoff = pw_grid%cutoff
207 IF (PRESENT(orthorhombic)) orthorhombic = pw_grid%orthorhombic
208 IF (PRESENT(gvectors)) gvectors => pw_grid%g
209 IF (PRESENT(gsquare)) gsquare => pw_grid%gsq
210
211 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 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 cpassert(pw_grid%ref_count > 0)
235
236 IF (PRESENT(grid_span)) THEN
237 pw_grid%grid_span = grid_span
238 END IF
239 IF (PRESENT(bounds) .AND. PRESENT(npts)) THEN
240 pw_grid%bounds = bounds
241 pw_grid%npts = npts
242 cpassert(all(npts == bounds(2, :) - bounds(1, :) + 1))
243 ELSE IF (PRESENT(bounds)) THEN
244 pw_grid%bounds = bounds
245 pw_grid%npts = bounds(2, :) - bounds(1, :) + 1
246 ELSE IF (PRESENT(npts)) THEN
247 pw_grid%npts = npts
248 pw_grid%bounds = pw_grid_bounds_from_n(npts)
249 END IF
250 IF (PRESENT(cutoff)) THEN
251 pw_grid%cutoff = cutoff
252 IF (PRESENT(spherical)) THEN
253 pw_grid%spherical = spherical
254 ELSE
255 pw_grid%spherical = .false.
256 END IF
257 END IF
258
259 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 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 CALL timeset(routinen, handle)
309
310 cell_deth = abs(det_3x3(cell_hmat))
311 IF (cell_deth < 1.0e-10_dp) THEN
312 CALL cp_abort(__location__, &
313 "An invalid set of cell vectors was specified. "// &
314 "The determinant det(h) is too small")
315 END IF
316 cell_h_inv = inv_3x3(cell_hmat)
317
318 IF (PRESENT(grid_span)) THEN
319 CALL set_pw_grid_info(pw_grid, grid_span=grid_span)
320 END IF
321
322 IF (PRESENT(spherical)) THEN
323 my_spherical = spherical
324 ELSE
325 my_spherical = .false.
326 END IF
327
328 IF (PRESENT(odd)) THEN
329 my_odd = odd
330 ELSE
331 my_odd = .false.
332 END IF
333
334 IF (PRESENT(fft_usage)) THEN
335 my_fft_usage = fft_usage
336 ELSE
337 my_fft_usage = .false.
338 END IF
339
340 IF (PRESENT(ncommensurate)) THEN
341 my_ncommensurate = ncommensurate
342 IF (PRESENT(icommensurate)) THEN
343 my_icommensurate = icommensurate
344 ELSE
345 my_icommensurate = min(1, ncommensurate)
346 END IF
347 ELSE
348 my_ncommensurate = 0
349 my_icommensurate = 1
350 END IF
351
352 IF (PRESENT(bounds)) THEN
353 IF (PRESENT(cutoff)) THEN
354 CALL set_pw_grid_info(pw_grid, bounds=bounds, cutoff=cutoff, &
355 spherical=my_spherical)
356 ELSE
357 n = bounds(2, :) - bounds(1, :) + 1
358 my_cutoff = pw_find_cutoff(n, cell_h_inv)
359 my_cutoff = 0.5_dp*my_cutoff*my_cutoff
360 CALL set_pw_grid_info(pw_grid, bounds=bounds, cutoff=my_cutoff, &
361 spherical=my_spherical)
362 END IF
363 ELSE IF (PRESENT(npts)) THEN
364 n = npts
365 IF (PRESENT(cutoff)) THEN
366 my_cutoff = cutoff
367 ELSE
368 my_cutoff = pw_find_cutoff(npts, cell_h_inv)
369 my_cutoff = 0.5_dp*my_cutoff*my_cutoff
370 END IF
371 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 ref_grid=ref_grid, n_orig=n)
376 END IF
377 CALL set_pw_grid_info(pw_grid, npts=n, cutoff=my_cutoff, &
378 spherical=my_spherical)
379 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 ref_grid=ref_grid)
384 CALL set_pw_grid_info(pw_grid, npts=n, cutoff=cutoff, &
385 spherical=my_spherical)
386 ELSE
387 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 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 CALL timestop(handle)
398
399 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 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 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: yz_mask
518 REAL(kind=dp) :: ecut
519
520!------------------------------------------------------------------------------
521
522 CALL timeset(routinen, handle)
523
524 cpassert(pw_grid%ref_count > 0)
525
526 ! set pointer to possible reference grid
527 IF (PRESENT(ref_grid)) THEN
528 pw_grid%reference = ref_grid%id_nr
529 END IF
530
531 IF (pw_grid%spherical) THEN
532 ecut = pw_grid%cutoff
533 ELSE
534 ecut = 1.e10_dp
535 END IF
536
537 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 ALLOCATE (yz_mask(n(2), n(3)))
544 CALL pw_grid_count(cell_h_inv, pw_grid, ecut, yz_mask)
545
546 ! Check if reference grid is compatible
547 IF (PRESENT(ref_grid)) THEN
548 cpassert(pw_grid%para%mode == ref_grid%para%mode)
549 cpassert(pw_grid%grid_span == ref_grid%grid_span)
550 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 rs_dims=rs_dims)
556
557 ! Allocate the grid fields
558 CALL pw_grid_allocate(pw_grid, pw_grid%ngpts_cut_local, &
559 pw_grid%bounds)
560
561 ! Fill in the grid structure
562 CALL pw_grid_assign(cell_h_inv, pw_grid, ecut)
563
564 ! Sort g vector wrt length (only local for each processor)
565 CALL pw_grid_sort(pw_grid, ref_grid)
566
567 CALL pw_grid_remap(pw_grid, yz_mask)
568
569 DEALLOCATE (yz_mask)
570
571 CALL cell2grid(cell_hmat, cell_h_inv, cell_deth, pw_grid)
572 !
573 ! Output: All the information of this grid type
574 !
575
576 IF (PRESENT(iounit)) THEN
577 CALL pw_grid_print(pw_grid, iounit)
578 END IF
579
580 CALL timestop(handle)
581
582 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 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 pw_grid%vol = abs(cell_deth)
599 pw_grid%dvol = pw_grid%vol/real(pw_grid%ngpts, kind=dp)
600 pw_grid%dr(1) = sqrt(sum(cell_hmat(:, 1)**2)) &
601 /real(pw_grid%npts(1), kind=dp)
602 pw_grid%dr(2) = sqrt(sum(cell_hmat(:, 2)**2)) &
603 /real(pw_grid%npts(2), kind=dp)
604 pw_grid%dr(3) = sqrt(sum(cell_hmat(:, 3)**2)) &
605 /real(pw_grid%npts(3), kind=dp)
606 pw_grid%dh(:, 1) = cell_hmat(:, 1)/real(pw_grid%npts(1), kind=dp)
607 pw_grid%dh(:, 2) = cell_hmat(:, 2)/real(pw_grid%npts(2), kind=dp)
608 pw_grid%dh(:, 3) = cell_hmat(:, 3)/real(pw_grid%npts(3), kind=dp)
609 pw_grid%dh_inv(1, :) = cell_h_inv(1, :)*real(pw_grid%npts(1), kind=dp)
610 pw_grid%dh_inv(2, :) = cell_h_inv(2, :)*real(pw_grid%npts(2), kind=dp)
611 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 (cell_hmat(3, 1) == 0.0_dp) .AND. (cell_hmat(3, 2) == 0.0_dp)) THEN
616 pw_grid%orthorhombic = .true.
617 ELSE
618 pw_grid%orthorhombic = .false.
619 END IF
620 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 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 IF (pw_grid%para%mode == pw_mode_local) THEN
643 IF (info > 0) THEN
644 WRITE (info, '(/,A,T71,I10)') &
645 " PW_GRID| Information for grid number ", pw_grid%id_nr
646 IF (pw_grid%reference > 0) THEN
647 WRITE (info, '(A,T71,I10)') &
648 " PW_GRID| Number of the reference grid ", pw_grid%reference
649 END IF
650 WRITE (info, '(" PW_GRID| Cutoff [a.u.]",T71,f10.1)') pw_grid%cutoff
651 IF (pw_grid%spherical) THEN
652 WRITE (info, '(A,T78,A)') " PW_GRID| spherical cutoff: ", "YES"
653 WRITE (info, '(A,T71,I10)') " PW_GRID| Grid points within cutoff", &
654 pw_grid%ngpts_cut
655 ELSE
656 WRITE (info, '(A,T78,A)') " PW_GRID| spherical cutoff: ", " NO"
657 END IF
658 DO i = 1, 3
659 WRITE (info, '(A,I3,T30,2I8,T62,A,T71,I10)') " PW_GRID| Bounds ", &
660 i, pw_grid%bounds(1, i), pw_grid%bounds(2, i), &
661 "Points:", pw_grid%npts(i)
662 END DO
663 WRITE (info, '(A,G12.4,T50,A,T67,F14.4)') &
664 " PW_GRID| Volume element (a.u.^3)", &
665 pw_grid%dvol, " Volume (a.u.^3) :", pw_grid%vol
666 IF (pw_grid%grid_span == halfspace) THEN
667 WRITE (info, '(A,T72,A)') " PW_GRID| Grid span", "HALFSPACE"
668 ELSE
669 WRITE (info, '(A,T72,A)') " PW_GRID| Grid span", "FULLSPACE"
670 END IF
671 END IF
672 ELSE
673
674 n(1) = pw_grid%ngpts_cut_local
675 n(2) = pw_grid%ngpts_local
676 CALL pw_grid%para%group%sum(n(1:2))
677 n(3) = sum(pw_grid%para%nyzray)
678 rv(:, 1) = real(n, kind=dp)/real(pw_grid%para%group_size, kind=dp)
679 n(1) = pw_grid%ngpts_cut_local
680 n(2) = pw_grid%ngpts_local
681 CALL pw_grid%para%group%max(n(1:2))
682 n(3) = maxval(pw_grid%para%nyzray)
683 rv(:, 2) = real(n, kind=dp)
684 n(1) = pw_grid%ngpts_cut_local
685 n(2) = pw_grid%ngpts_local
686 CALL pw_grid%para%group%min(n(1:2))
687 n(3) = minval(pw_grid%para%nyzray)
688 rv(:, 3) = real(n, kind=dp)
689
690 IF (pw_grid%para%group_head .AND. info > 0) THEN
691 WRITE (info, '(/,A,T71,I10)') &
692 " PW_GRID| Information for grid number ", pw_grid%id_nr
693 IF (pw_grid%reference > 0) THEN
694 WRITE (info, '(A,T71,I10)') &
695 " PW_GRID| Number of the reference grid ", pw_grid%reference
696 END IF
697 WRITE (info, '(A,T60,I10,A)') &
698 " PW_GRID| Grid distributed over ", pw_grid%para%group_size, &
699 " processors"
700 WRITE (info, '(A,T71,2I5)') &
701 " PW_GRID| Real space group dimensions ", pw_grid%para%rs_dims
702 IF (pw_grid%para%blocked) THEN
703 WRITE (info, '(A,T78,A)') " PW_GRID| the grid is blocked: ", "YES"
704 ELSE
705 WRITE (info, '(A,T78,A)') " PW_GRID| the grid is blocked: ", " NO"
706 END IF
707 WRITE (info, '(" PW_GRID| Cutoff [a.u.]",T71,f10.1)') pw_grid%cutoff
708 IF (pw_grid%spherical) THEN
709 WRITE (info, '(A,T78,A)') " PW_GRID| spherical cutoff: ", "YES"
710 WRITE (info, '(A,T71,I10)') " PW_GRID| Grid points within cutoff", &
711 pw_grid%ngpts_cut
712 ELSE
713 WRITE (info, '(A,T78,A)') " PW_GRID| spherical cutoff: ", " NO"
714 END IF
715 DO i = 1, 3
716 WRITE (info, '(A,I3,T30,2I8,T62,A,T71,I10)') " PW_GRID| Bounds ", &
717 i, pw_grid%bounds(1, i), pw_grid%bounds(2, i), &
718 "Points:", pw_grid%npts(i)
719 END DO
720 WRITE (info, '(A,G12.4,T50,A,T67,F14.4)') &
721 " PW_GRID| Volume element (a.u.^3)", &
722 pw_grid%dvol, " Volume (a.u.^3) :", pw_grid%vol
723 IF (pw_grid%grid_span == halfspace) THEN
724 WRITE (info, '(A,T72,A)') " PW_GRID| Grid span", "HALFSPACE"
725 ELSE
726 WRITE (info, '(A,T72,A)') " PW_GRID| Grid span", "FULLSPACE"
727 END IF
728 WRITE (info, '(A,T48,A)') " PW_GRID| Distribution", &
729 " Average Max Min"
730 WRITE (info, '(A,T45,F12.1,2I12)') " PW_GRID| G-Vectors", &
731 rv(1, 1), nint(rv(1, 2)), nint(rv(1, 3))
732 WRITE (info, '(A,T45,F12.1,2I12)') " PW_GRID| G-Rays ", &
733 rv(3, 1), nint(rv(3, 2)), nint(rv(3, 3))
734 WRITE (info, '(A,T45,F12.1,2I12)') " PW_GRID| Real Space Points", &
735 rv(2, 1), nint(rv(2, 2)), nint(rv(2, 3))
736 END IF ! group head
737 END IF ! local
738
739 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: pemap
771 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: yz_index
772 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 CALL timeset(routinen, handle)
780
781 lby = pw_grid%bounds(1, 2)
782 lbz = pw_grid%bounds(1, 3)
783
784 pw_grid%ngpts = product(int(pw_grid%npts, kind=int_8))
785 cpassert(all(pw_grid%para%rs_dims == 0))
786 IF (PRESENT(rs_dims)) THEN
787 pw_grid%para%rs_dims = rs_dims
788 END IF
789
790 IF (PRESENT(blocked)) THEN
791 blocked_local = blocked
792 ELSE
793 blocked_local = do_pw_grid_blocked_free
794 END IF
795
796 pw_grid%para%blocked = .false.
797
798 IF (pw_grid%para%mode == pw_mode_local) THEN
799
800 pw_grid%para%ray_distribution = .false.
801
802 pw_grid%bounds_local = pw_grid%bounds
803 pw_grid%npts_local = pw_grid%npts
804 cpassert(pw_grid%ngpts_cut < huge(pw_grid%ngpts_cut_local))
805 pw_grid%ngpts_cut_local = int(pw_grid%ngpts_cut)
806 cpassert(pw_grid%ngpts < huge(pw_grid%ngpts_local))
807 pw_grid%ngpts_local = int(pw_grid%ngpts)
808 pw_grid%para%rs_dims = 1
809 mp_comm_old = mp_comm_self
810 CALL pw_grid%para%rs_group%create(mp_comm_old, 2, &
811 pw_grid%para%rs_dims)
812 pw_grid%para%rs_dims = pw_grid%para%rs_group%num_pe_cart
813 pw_grid%para%rs_pos = pw_grid%para%rs_group%mepos_cart
814 CALL pw_grid%para%rs_group%rank_cart(pw_grid%para%rs_pos, pw_grid%para%rs_mpo)
815 pw_grid%para%group_size = 1
816 ALLOCATE (pw_grid%para%bo(2, 3, 0:0, 3))
817 DO i = 1, 3
818 pw_grid%para%bo(1, 1:3, 0, i) = 1
819 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 nx = pw_grid%npts(1)
826 ny = pw_grid%npts(2)
827 nz = pw_grid%npts(3)
828 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 IF (any(pw_grid%para%rs_dims <= 0)) THEN
840 IF (all(pw_grid%para%rs_dims <= 0)) THEN
841 pw_grid%para%rs_dims = (/0, 0/)
842 ELSE
843 IF (pw_grid%para%rs_dims(1) > 0) THEN
844 pw_grid%para%rs_dims(2) = np/pw_grid%para%rs_dims(1)
845 ELSE
846 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 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 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 IF (all(pw_grid%para%rs_dims == (/0, 0/))) THEN
857 ! only small grids have a chance to be 2d distributed
858 IF (nx < np) THEN
859 ! gives the most square looking distribution
860 CALL mp_dims_create(np, pw_grid%para%rs_dims)
861 ! we tend to like the first index being smaller than the second
862 IF (pw_grid%para%rs_dims(1) > pw_grid%para%rs_dims(2)) THEN
863 itmp = pw_grid%para%rs_dims(1)
864 pw_grid%para%rs_dims(1) = pw_grid%para%rs_dims(2)
865 pw_grid%para%rs_dims(2) = itmp
866 END IF
867 ! but should avoid having the first index 1 in all cases
868 IF (pw_grid%para%rs_dims(1) == 1) THEN
869 itmp = pw_grid%para%rs_dims(1)
870 pw_grid%para%rs_dims(1) = pw_grid%para%rs_dims(2)
871 pw_grid%para%rs_dims(2) = itmp
872 END IF
873 ELSE
874 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 blocking = .false.
882 CASE (do_pw_grid_blocked_true)
883 blocking = .true.
884 CASE (do_pw_grid_blocked_free)
885 IF (all(pw_grid%para%rs_dims == (/np, 1/))) THEN
886 blocking = .false.
887 ELSE
888 blocking = .true.
889 END IF
890 CASE DEFAULT
891 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 pw_grid%para%rs_dims)
897 pw_grid%para%rs_dims = pw_grid%para%rs_group%num_pe_cart
898 pw_grid%para%rs_pos = pw_grid%para%rs_group%mepos_cart
899 CALL pw_grid%para%rs_group%rank_cart(pw_grid%para%rs_pos, pw_grid%para%rs_mpo)
900
901 IF (PRESENT(bounds_local)) THEN
902 pw_grid%bounds_local = bounds_local
903 ELSE
904 lo = get_limit(nx, pw_grid%para%rs_dims(1), &
905 pw_grid%para%rs_pos(1))
906 pw_grid%bounds_local(:, 1) = lo + pw_grid%bounds(1, 1) - 1
907 lo = get_limit(ny, pw_grid%para%rs_dims(2), &
908 pw_grid%para%rs_pos(2))
909 pw_grid%bounds_local(:, 2) = lo + pw_grid%bounds(1, 2) - 1
910 pw_grid%bounds_local(:, 3) = pw_grid%bounds(:, 3)
911 END IF
912
913 pw_grid%npts_local(:) = pw_grid%bounds_local(2, :) &
914 - pw_grid%bounds_local(1, :) + 1
915
916 !..the third distribution is needed for the second step in the FFT
917 ALLOCATE (pw_grid%para%bo(2, 3, 0:np - 1, 3))
918 rsd = pw_grid%para%rs_dims
919
920 IF (PRESENT(bounds_local)) THEN
921 ! axis_dist tells what portion of 1 .. nx , 1 .. ny , 1 .. nz are in the current process
922 DO i = 1, 3
923 axis_dist(:, i) = bounds_local(:, i) - pw_grid%bounds(1, i) + 1
924 END DO
925 ALLOCATE (axis_dist_all(2, 3, np))
926 CALL pw_grid%para%rs_group%allgather(axis_dist, axis_dist_all)
927 DO ip = 0, np - 1
928 CALL pw_grid%para%rs_group%coords(ip, coor)
929 ! distribution xyZ
930 pw_grid%para%bo(1:2, 1, ip, 1) = axis_dist_all(1:2, 1, ip + 1)
931 pw_grid%para%bo(1:2, 2, ip, 1) = axis_dist_all(1:2, 2, ip + 1)
932 pw_grid%para%bo(1, 3, ip, 1) = 1
933 pw_grid%para%bo(2, 3, ip, 1) = nz
934 ! distribution xYz
935 pw_grid%para%bo(1:2, 1, ip, 2) = axis_dist_all(1:2, 1, ip + 1)
936 pw_grid%para%bo(1, 2, ip, 2) = 1
937 pw_grid%para%bo(2, 2, ip, 2) = ny
938 pw_grid%para%bo(1:2, 3, ip, 2) = get_limit(nz, rsd(2), coor(2))
939 ! distribution Xyz
940 pw_grid%para%bo(1, 1, ip, 3) = 1
941 pw_grid%para%bo(2, 1, ip, 3) = nx
942 pw_grid%para%bo(1:2, 2, ip, 3) = get_limit(ny, rsd(1), coor(1))
943 pw_grid%para%bo(1:2, 3, ip, 3) = get_limit(nz, rsd(2), coor(2))
944 END DO
945 DEALLOCATE (axis_dist_all)
946 ELSE
947 DO ip = 0, np - 1
948 CALL pw_grid%para%rs_group%coords(ip, coor)
949 ! distribution xyZ
950 pw_grid%para%bo(1:2, 1, ip, 1) = get_limit(nx, rsd(1), coor(1))
951 pw_grid%para%bo(1:2, 2, ip, 1) = get_limit(ny, rsd(2), coor(2))
952 pw_grid%para%bo(1, 3, ip, 1) = 1
953 pw_grid%para%bo(2, 3, ip, 1) = nz
954 ! distribution xYz
955 pw_grid%para%bo(1:2, 1, ip, 2) = get_limit(nx, rsd(1), coor(1))
956 pw_grid%para%bo(1, 2, ip, 2) = 1
957 pw_grid%para%bo(2, 2, ip, 2) = ny
958 pw_grid%para%bo(1:2, 3, ip, 2) = get_limit(nz, rsd(2), coor(2))
959 ! distribution Xyz
960 pw_grid%para%bo(1, 1, ip, 3) = 1
961 pw_grid%para%bo(2, 1, ip, 3) = nx
962 pw_grid%para%bo(1:2, 2, ip, 3) = get_limit(ny, rsd(1), coor(1))
963 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 pw_grid%ngpts_cut_local = 0
968
969 ALLOCATE (pw_grid%para%nyzray(0:np - 1))
970
971 ALLOCATE (pw_grid%para%yzq(ny, nz))
972
973 IF (pw_grid%spherical .OR. pw_grid%grid_span == halfspace &
974 .OR. .NOT. blocking) THEN
975
976 pw_grid%para%ray_distribution = .true.
977
978 pw_grid%para%yzq = -1
979 IF (PRESENT(ref_grid)) THEN
980 ! tag all vectors from the reference grid
981 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 i1 = SIZE(yz_mask, 1)
989 i2 = SIZE(yz_mask, 2)
990 ALLOCATE (yz_index(2, i1*i2))
991 CALL order_mask(yz_mask, yz_index)
992 DO i = 1, i1*i2
993 lo(1) = yz_index(1, i)
994 lo(2) = yz_index(2, i)
995 IF (lo(1)*lo(2) == 0) cycle
996 gmax = yz_mask(lo(1), lo(2))
997 IF (gmax == 0) cycle
998 yz_mask(lo(1), lo(2)) = 0
999 ip = mod(i - 1, 2*np)
1000 IF (ip > np - 1) ip = 2*np - ip - 1
1001 IF (ip == pw_grid%para%my_pos) THEN
1002 pw_grid%ngpts_cut_local = pw_grid%ngpts_cut_local + gmax
1003 END IF
1004 pw_grid%para%yzq(lo(1), lo(2)) = ip
1005 IF (pw_grid%grid_span == halfspace) THEN
1006 m = -lo(1) - 2*lby + 2
1007 n = -lo(2) - 2*lbz + 2
1008 pw_grid%para%yzq(m, n) = ip
1009 yz_mask(m, n) = 0
1010 END IF
1011 END DO
1012
1013 DEALLOCATE (yz_index)
1014
1015 ! Count the total number of rays on each processor
1016 pw_grid%para%nyzray = 0
1017 DO i = 1, nz
1018 DO j = 1, ny
1019 ip = pw_grid%para%yzq(j, i)
1020 IF (ip >= 0) pw_grid%para%nyzray(ip) = &
1021 pw_grid%para%nyzray(ip) + 1
1022 END DO
1023 END DO
1024
1025 ! Allocate mapping array (y:z, nray, nproc)
1026 ns = maxval(pw_grid%para%nyzray(0:np - 1))
1027 ALLOCATE (pw_grid%para%yzp(2, ns, 0:np - 1))
1028
1029 ! Fill mapping array, recalculate nyzray for convenience
1030 pw_grid%para%nyzray = 0
1031 DO i = 1, nz
1032 DO j = 1, ny
1033 ip = pw_grid%para%yzq(j, i)
1034 IF (ip >= 0) THEN
1035 pw_grid%para%nyzray(ip) = &
1036 pw_grid%para%nyzray(ip) + 1
1037 ns = pw_grid%para%nyzray(ip)
1038 pw_grid%para%yzp(1, ns, ip) = j
1039 pw_grid%para%yzp(2, ns, ip) = i
1040 IF (ip == pw_grid%para%my_pos) THEN
1041 pw_grid%para%yzq(j, i) = ns
1042 ELSE
1043 pw_grid%para%yzq(j, i) = -1
1044 END IF
1045 ELSE
1046 pw_grid%para%yzq(j, i) = -2
1047 END IF
1048 END DO
1049 END DO
1050
1051 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 pw_grid%para%blocked = .true.
1059 pw_grid%para%ray_distribution = .false.
1060
1061 DO ip = 0, np - 1
1062 m = pw_grid%para%bo(2, 2, ip, 3) - &
1063 pw_grid%para%bo(1, 2, ip, 3) + 1
1064 n = pw_grid%para%bo(2, 3, ip, 3) - &
1065 pw_grid%para%bo(1, 3, ip, 3) + 1
1066 pw_grid%para%nyzray(ip) = n*m
1067 END DO
1068
1069 ipl = pw_grid%para%rs_mpo
1070 l = pw_grid%para%bo(2, 1, ipl, 3) - &
1071 pw_grid%para%bo(1, 1, ipl, 3) + 1
1072 m = pw_grid%para%bo(2, 2, ipl, 3) - &
1073 pw_grid%para%bo(1, 2, ipl, 3) + 1
1074 n = pw_grid%para%bo(2, 3, ipl, 3) - &
1075 pw_grid%para%bo(1, 3, ipl, 3) + 1
1076 pw_grid%ngpts_cut_local = l*m*n
1077 pw_grid%ngpts_local = pw_grid%ngpts_cut_local
1078
1079 pw_grid%para%yzq = 0
1080 ny = pw_grid%para%bo(2, 2, ipl, 3) - &
1081 pw_grid%para%bo(1, 2, ipl, 3) + 1
1082 DO n = pw_grid%para%bo(1, 3, ipl, 3), &
1083 pw_grid%para%bo(2, 3, ipl, 3)
1084 i = n - pw_grid%para%bo(1, 3, ipl, 3)
1085 DO m = pw_grid%para%bo(1, 2, ipl, 3), &
1086 pw_grid%para%bo(2, 2, ipl, 3)
1087 j = m - pw_grid%para%bo(1, 2, ipl, 3) + 1
1088 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 ns = maxval(pw_grid%para%nyzray(0:np - 1))
1094 ALLOCATE (pw_grid%para%yzp(2, ns, 0:np - 1))
1095 pw_grid%para%yzp = 0
1096
1097 ALLOCATE (pemap(0:np - 1))
1098 pemap = 0
1099 pemap(pw_grid%para%my_pos) = pw_grid%para%rs_mpo
1100 CALL pw_grid%para%group%sum(pemap)
1101
1102 DO ip = 0, np - 1
1103 ipp = pemap(ip)
1104 ns = 0
1105 DO n = pw_grid%para%bo(1, 3, ipp, 3), &
1106 pw_grid%para%bo(2, 3, ipp, 3)
1107 i = n - pw_grid%bounds(1, 3) + 1
1108 DO m = pw_grid%para%bo(1, 2, ipp, 3), &
1109 pw_grid%para%bo(2, 2, ipp, 3)
1110 j = m - pw_grid%bounds(1, 2) + 1
1111 ns = ns + 1
1112 pw_grid%para%yzp(1, ns, ip) = j
1113 pw_grid%para%yzp(2, ns, ip) = i
1114 END DO
1115 END DO
1116 cpassert(ns == pw_grid%para%nyzray(ip))
1117 END DO
1118
1119 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 IF (pw_grid%para%mode .EQ. pw_mode_distributed) THEN
1128 ALLOCATE (pw_grid%para%pos_of_x(pw_grid%bounds(1, 1):pw_grid%bounds(2, 1)))
1129 pw_grid%para%pos_of_x = 0
1130 pw_grid%para%pos_of_x(pw_grid%bounds_local(1, 1):pw_grid%bounds_local(2, 1)) = pw_grid%para%my_pos
1131 CALL pw_grid%para%group%sum(pw_grid%para%pos_of_x)
1132 ELSE
1133 ! this should not be needed
1134 ALLOCATE (pw_grid%para%pos_of_x(pw_grid%bounds(1, 1):pw_grid%bounds(2, 1)))
1135 pw_grid%para%pos_of_x = 0
1136 END IF
1137
1138 CALL timestop(handle)
1139
1140 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 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 ny = ref_grid%npts(2)
1160 nz = ref_grid%npts(3)
1161 lby = pw_grid%bounds(1, 2)
1162 lbz = pw_grid%bounds(1, 3)
1163 uby = pw_grid%bounds(2, 2)
1164 ubz = pw_grid%bounds(2, 3)
1165 my = SIZE(yz_mask, 1)
1166 mz = SIZE(yz_mask, 2)
1167
1168 ! loop over all processors and all g vectors yz lines on this processor
1169 DO ip = 0, ref_grid%para%group_size - 1
1170 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 y = ref_grid%para%yzp(1, ig, ip) - 1
1174 IF (y >= ny/2) y = y - ny
1175 z = ref_grid%para%yzp(2, ig, ip) - 1
1176 IF (z >= nz/2) z = z - nz
1177 ! check if this is inside the realm of the new grid
1178 IF (y < lby .OR. y > uby .OR. z < lbz .OR. z > ubz) cycle
1179 ! go to shifted coordinates
1180 y = y - lby + 1
1181 z = z - lbz + 1
1182 ! this tag is outside the cutoff range of the new grid
1183 IF (pw_grid%grid_span == halfspace) THEN
1184 yp = -y - 2*lby + 2
1185 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 IF (yp < 1 .OR. yp > my .OR. zp < 1 .OR. zp > mz) cycle
1189 gmax = max(yz_mask(y, z), yz_mask(yp, zp))
1190 IF (gmax == 0) cycle
1191 yz_mask(y, z) = 0
1192 yz_mask(yp, zp) = 0
1193 pw_grid%para%yzq(y, z) = ip
1194 pw_grid%para%yzq(yp, zp) = ip
1195 ELSE
1196 gmax = yz_mask(y, z)
1197 IF (gmax == 0) cycle
1198 yz_mask(y, z) = 0
1199 pw_grid%para%yzq(y, z) = ip
1200 END IF
1201 IF (ip == pw_grid%para%my_pos) THEN
1202 pw_grid%ngpts_cut_local = pw_grid%ngpts_cut_local + gmax
1203 END IF
1204 END DO
1205 END DO
1206
1207 END SUBROUTINE pre_tag
1208
1209! **************************************************************************************************
1210!> \brief ...
1211!> \param yz_mask ...
1212!> \param yz_index ...
1213! **************************************************************************************************
1214 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 i1 = SIZE(yz_mask, 1)
1228 i2 = SIZE(yz_mask, 2)
1229 yz_index = 0
1230
1231 icount = 1
1232 ic = i1/2
1233 jc = i2/2
1234 ii = ic
1235 jj = jc
1236 IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1237 IF (yz_mask(ii, jj) /= 0) THEN
1238 yz_index(1, icount) = ii
1239 yz_index(2, icount) = jj
1240 icount = icount + 1
1241 END IF
1242 END IF
1243 DO im = 1, max(ic + 1, jc + 1)
1244 ii = ic - im
1245 DO jj = jc - im, jc + im
1246 IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1247 IF (yz_mask(ii, jj) /= 0) THEN
1248 yz_index(1, icount) = ii
1249 yz_index(2, icount) = jj
1250 icount = icount + 1
1251 END IF
1252 END IF
1253 END DO
1254 ii = ic + im
1255 DO jj = jc - im, jc + im
1256 IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1257 IF (yz_mask(ii, jj) /= 0) THEN
1258 yz_index(1, icount) = ii
1259 yz_index(2, icount) = jj
1260 icount = icount + 1
1261 END IF
1262 END IF
1263 END DO
1264 jj = jc - im
1265 DO ii = ic - im + 1, ic + im - 1
1266 IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1267 IF (yz_mask(ii, jj) /= 0) THEN
1268 yz_index(1, icount) = ii
1269 yz_index(2, icount) = jj
1270 icount = icount + 1
1271 END IF
1272 END IF
1273 END DO
1274 jj = jc + im
1275 DO ii = ic - im + 1, ic + im - 1
1276 IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN
1277 IF (yz_mask(ii, jj) /= 0) THEN
1278 yz_index(1, icount) = ii
1279 yz_index(2, icount) = jj
1280 icount = icount + 1
1281 END IF
1282 END IF
1283 END DO
1284 END DO
1285
1286 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 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 + 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 + 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 + real(n, dp)*h_inv(3, 3)
1316
1317 ! enforce strict zero-ness in this case (compiler optimization)
1318 IF (l == 0 .AND. m == 0 .AND. n == 0) THEN
1319 length_x = 0
1320 length_y = 0
1321 length_z = 0
1322 END IF
1323
1324 length_x = length_x*twopi
1325 length_y = length_y*twopi
1326 length_z = length_z*twopi
1327
1328 length = length_x**2 + length_y**2 + length_z**2
1329
1330 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 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 IF (pw_grid%grid_span == halfspace) THEN
1357 n_upperlimit = 0
1358 ELSE IF (pw_grid%grid_span == fullspace) THEN
1359 n_upperlimit = bounds(2, 3)
1360 ELSE
1361 cpabort("No type set for the grid")
1362 END IF
1363
1364 ! finds valid g-points within grid
1365 gpt = 0
1366 IF (pw_grid%para%mode == pw_mode_local) THEN
1367 nlim(1) = bounds(1, 3)
1368 nlim(2) = n_upperlimit
1369 ELSE IF (pw_grid%para%mode == pw_mode_distributed) THEN
1370 n = n_upperlimit - bounds(1, 3) + 1
1371 nlim = get_limit(n, pw_grid%para%group_size, pw_grid%para%my_pos)
1372 nlim = nlim + bounds(1, 3) - 1
1373 ELSE
1374 cpabort("para % mode not specified")
1375 END IF
1376
1377 yz_mask = 0
1378 DO n = nlim(1), nlim(2)
1379 nn = n - bounds(1, 3) + 1
1380 DO m = bounds(1, 2), bounds(2, 2)
1381 mm = m - bounds(1, 2) + 1
1382 DO l = bounds(1, 1), bounds(2, 1)
1383 IF (pw_grid%grid_span == halfspace .AND. n == 0) THEN
1384 IF ((m == 0 .AND. l > 0) .OR. (m > 0)) cycle
1385 END IF
1386
1387 CALL pw_vec_length(h_inv, length_x, length_y, length_z, length, l, m, n)
1388
1389 IF (0.5_dp*length <= cutoff) THEN
1390 gpt = gpt + 1
1391 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 IF (pw_grid%para%mode == pw_mode_distributed) THEN
1401 CALL pw_grid%para%group%sum(gpt)
1402 CALL pw_grid%para%group%sum(yz_mask)
1403 END IF
1404 pw_grid%ngpts_cut = gpt
1405
1406 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 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 CALL timeset(routinen, handle)
1433
1434 bounds = pw_grid%bounds
1435 lby = pw_grid%bounds(1, 2)
1436 lbz = pw_grid%bounds(1, 3)
1437
1438 IF (pw_grid%grid_span == halfspace) THEN
1439 n_upperlimit = 0
1440 ELSE IF (pw_grid%grid_span == fullspace) THEN
1441 n_upperlimit = bounds(2, 3)
1442 ELSE
1443 cpabort("No type set for the grid")
1444 END IF
1445
1446 ! finds valid g-points within grid
1447 IF (pw_grid%para%mode == pw_mode_local) THEN
1448 gpt = 0
1449 DO n = bounds(1, 3), n_upperlimit
1450 DO m = bounds(1, 2), bounds(2, 2)
1451 DO l = bounds(1, 1), bounds(2, 1)
1452 IF (pw_grid%grid_span == halfspace .AND. n == 0) THEN
1453 IF ((m == 0 .AND. l > 0) .OR. (m > 0)) cycle
1454 END IF
1455
1456 CALL pw_vec_length(h_inv, length_x, length_y, length_z, length, l, m, n)
1457
1458 IF (0.5_dp*length <= cutoff) THEN
1459 gpt = gpt + 1
1460 pw_grid%g(1, gpt) = length_x
1461 pw_grid%g(2, gpt) = length_y
1462 pw_grid%g(3, gpt) = length_z
1463 pw_grid%gsq(gpt) = length
1464 pw_grid%g_hat(1, gpt) = l
1465 pw_grid%g_hat(2, gpt) = m
1466 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 IF (pw_grid%para%ray_distribution) THEN
1476
1477 gpt = 0
1478 ip = pw_grid%para%my_pos
1479 DO i = 1, pw_grid%para%nyzray(ip)
1480 n = pw_grid%para%yzp(2, i, ip) + lbz - 1
1481 m = pw_grid%para%yzp(1, i, ip) + lby - 1
1482 IF (n > n_upperlimit) cycle
1483 DO l = bounds(1, 1), bounds(2, 1)
1484 IF (pw_grid%grid_span == halfspace .AND. n == 0) THEN
1485 IF ((m == 0 .AND. l > 0) .OR. (m > 0)) cycle
1486 END IF
1487
1488 CALL pw_vec_length(h_inv, length_x, length_y, length_z, length, l, m, n)
1489
1490 IF (0.5_dp*length <= cutoff) THEN
1491 gpt = gpt + 1
1492 pw_grid%g(1, gpt) = length_x
1493 pw_grid%g(2, gpt) = length_y
1494 pw_grid%g(3, gpt) = length_z
1495 pw_grid%gsq(gpt) = length
1496 pw_grid%g_hat(1, gpt) = l
1497 pw_grid%g_hat(2, gpt) = m
1498 pw_grid%g_hat(3, gpt) = n
1499 END IF
1500
1501 END DO
1502 END DO
1503
1504 ELSE
1505
1506 bol = pw_grid%para%bo(:, :, pw_grid%para%rs_mpo, 3)
1507 gpt = 0
1508 DO n = bounds(1, 3), bounds(2, 3)
1509 IF (n < 0) THEN
1510 nn = n + pw_grid%npts(3) + 1
1511 ELSE
1512 nn = n + 1
1513 END IF
1514 IF (nn < bol(1, 3) .OR. nn > bol(2, 3)) cycle
1515 DO m = bounds(1, 2), bounds(2, 2)
1516 IF (m < 0) THEN
1517 mm = m + pw_grid%npts(2) + 1
1518 ELSE
1519 mm = m + 1
1520 END IF
1521 IF (mm < bol(1, 2) .OR. mm > bol(2, 2)) cycle
1522 DO l = bounds(1, 1), bounds(2, 1)
1523 IF (l < 0) THEN
1524 ll = l + pw_grid%npts(1) + 1
1525 ELSE
1526 ll = l + 1
1527 END IF
1528 IF (ll < bol(1, 1) .OR. ll > bol(2, 1)) cycle
1529
1530 CALL pw_vec_length(h_inv, length_x, length_y, length_z, length, l, m, n)
1531
1532 gpt = gpt + 1
1533 pw_grid%g(1, gpt) = length_x
1534 pw_grid%g(2, gpt) = length_y
1535 pw_grid%g(3, gpt) = length_z
1536 pw_grid%gsq(gpt) = length
1537 pw_grid%g_hat(1, gpt) = l
1538 pw_grid%g_hat(2, gpt) = m
1539 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 cpassert(pw_grid%ngpts_cut_local == gpt)
1551 IF (pw_grid%para%mode == pw_mode_distributed) THEN
1552 gpt_global = gpt
1553 CALL pw_grid%para%group%sum(gpt_global)
1554 cpassert(pw_grid%ngpts_cut == gpt_global)
1555 END IF
1556
1557 pw_grid%have_g0 = .false.
1558 pw_grid%first_gne0 = 1
1559 DO gpt = 1, pw_grid%ngpts_cut_local
1560 IF (all(pw_grid%g_hat(:, gpt) == 0)) THEN
1561 pw_grid%have_g0 = .true.
1562 pw_grid%first_gne0 = 2
1563 EXIT
1564 END IF
1565 END DO
1566
1567 CALL pw_grid_set_maps(pw_grid%grid_span, pw_grid%g_hat, &
1568 pw_grid%mapl, pw_grid%mapm, pw_grid%mapn, pw_grid%npts)
1569
1570 CALL timestop(handle)
1571
1572 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 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 ng = SIZE(g_hat, 2)
1599
1600 DO gpt = 1, ng
1601 l = g_hat(1, gpt)
1602 m = g_hat(2, gpt)
1603 n = g_hat(3, gpt)
1604 IF (l < 0) THEN
1605 mapl%pos(l) = l + npts(1)
1606 ELSE
1607 mapl%pos(l) = l
1608 END IF
1609 IF (m < 0) THEN
1610 mapm%pos(m) = m + npts(2)
1611 ELSE
1612 mapm%pos(m) = m
1613 END IF
1614 IF (n < 0) THEN
1615 mapn%pos(n) = n + npts(3)
1616 ELSE
1617 mapn%pos(n) = n
1618 END IF
1619
1620 ! Generating the maps to the full 3-d space
1621
1622 IF (grid_span == halfspace) THEN
1623
1624 IF (l <= 0) THEN
1625 mapl%neg(l) = -l
1626 ELSE
1627 mapl%neg(l) = npts(1) - l
1628 END IF
1629 IF (m <= 0) THEN
1630 mapm%neg(m) = -m
1631 ELSE
1632 mapm%neg(m) = npts(2) - m
1633 END IF
1634 IF (n <= 0) THEN
1635 mapn%neg(n) = -n
1636 ELSE
1637 mapn%neg(n) = npts(3) - n
1638 END IF
1639
1640 END IF
1641
1642 END DO
1643
1644 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 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 CALL timeset(routinen, handle)
1677
1678 ALLOCATE (pw_grid%g(3, ng))
1679 ALLOCATE (pw_grid%gsq(ng))
1680 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 ALLOCATE (pw_grid%g_hatmap(1, 1))
1693#endif
1694
1695 IF (pw_grid%para%mode == pw_mode_distributed) THEN
1696 ALLOCATE (pw_grid%grays(pw_grid%npts(1), &
1697 pw_grid%para%nyzray(pw_grid%para%my_pos)))
1698 END IF
1699
1700 ALLOCATE (pw_grid%mapl%pos(bounds(1, 1):bounds(2, 1)))
1701 ALLOCATE (pw_grid%mapl%neg(bounds(1, 1):bounds(2, 1)))
1702 ALLOCATE (pw_grid%mapm%pos(bounds(1, 2):bounds(2, 2)))
1703 ALLOCATE (pw_grid%mapm%neg(bounds(1, 2):bounds(2, 2)))
1704 ALLOCATE (pw_grid%mapn%pos(bounds(1, 3):bounds(2, 3)))
1705 ALLOCATE (pw_grid%mapn%neg(bounds(1, 3):bounds(2, 3)))
1706
1707 CALL timestop(handle)
1708
1709 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 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 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: int_tmp
1741 LOGICAL :: g_found
1742 REAL(kind=dp) :: gig, gigr
1743 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: real_tmp
1744
1745 CALL timeset(routinen, handle)
1746
1747 ng = SIZE(pw_grid%gsq)
1748 ALLOCATE (idx(ng))
1749
1750 ! grids are (locally) ordered by length of G-vectors
1751 CALL sort(pw_grid%gsq, ng, idx)
1752 ! within shells order wrt x,y,z
1753 CALL sort_shells(pw_grid%gsq, pw_grid%g_hat, idx)
1754
1755 ALLOCATE (real_tmp(3, ng))
1756 DO i = 1, ng
1757 real_tmp(1, i) = pw_grid%g(1, idx(i))
1758 real_tmp(2, i) = pw_grid%g(2, idx(i))
1759 real_tmp(3, i) = pw_grid%g(3, idx(i))
1760 END DO
1761 DO i = 1, ng
1762 pw_grid%g(1, i) = real_tmp(1, i)
1763 pw_grid%g(2, i) = real_tmp(2, i)
1764 pw_grid%g(3, i) = real_tmp(3, i)
1765 END DO
1766 DEALLOCATE (real_tmp)
1767
1768 ALLOCATE (int_tmp(3, ng))
1769 DO i = 1, ng
1770 int_tmp(1, i) = pw_grid%g_hat(1, idx(i))
1771 int_tmp(2, i) = pw_grid%g_hat(2, idx(i))
1772 int_tmp(3, i) = pw_grid%g_hat(3, idx(i))
1773 END DO
1774 DO i = 1, ng
1775 pw_grid%g_hat(1, i) = int_tmp(1, i)
1776 pw_grid%g_hat(2, i) = int_tmp(2, i)
1777 pw_grid%g_hat(3, i) = int_tmp(3, i)
1778 END DO
1779 DEALLOCATE (int_tmp)
1780
1781 DEALLOCATE (idx)
1782
1783 ! check if ordering is compatible to reference grid
1784 IF (PRESENT(ref_grid)) THEN
1785 ngr = SIZE(ref_grid%gsq)
1786 ngr = min(ng, ngr)
1787 IF (pw_grid%spherical) THEN
1788 IF (.NOT. all(pw_grid%g_hat(1:3, 1:ngr) &
1789 == ref_grid%g_hat(1:3, 1:ngr))) THEN
1790 cpabort("G space sorting not compatible")
1791 END IF
1792 ELSE
1793 ALLOCATE (pw_grid%gidx(1:ngr))
1794 pw_grid%gidx = 0
1795 ! first try as many trivial associations as possible
1796 it = 0
1797 DO ig = 1, ngr
1798 IF (.NOT. all(pw_grid%g_hat(1:3, ig) &
1799 == ref_grid%g_hat(1:3, ig))) EXIT
1800 pw_grid%gidx(ig) = ig
1801 it = ig
1802 END DO
1803 ! now for the ones that could not be done
1804 IF (ng == ngr) THEN
1805 ! for the case pw_grid <= ref_grid
1806 is = it
1807 DO ig = it + 1, ngr
1808 gig = pw_grid%gsq(ig)
1809 gigr = max(1._dp, gig)
1810 g_found = .false.
1811 DO ih = is + 1, SIZE(ref_grid%gsq)
1812 IF (abs(gig - ref_grid%gsq(ih))/gigr > 1.e-12_dp) cycle
1813 g_found = .true.
1814 EXIT
1815 END DO
1816 IF (.NOT. g_found) THEN
1817 WRITE (*, "(A,I10,F20.10)") "G-vector", ig, pw_grid%gsq(ig)
1818 cpabort("G vector not found")
1819 END IF
1820 ip = ih - 1
1821 DO ih = ip + 1, SIZE(ref_grid%gsq)
1822 IF (abs(gig - ref_grid%gsq(ih))/gigr > 1.e-12_dp) cycle
1823 IF (pw_grid%g_hat(1, ig) /= ref_grid%g_hat(1, ih)) cycle
1824 IF (pw_grid%g_hat(2, ig) /= ref_grid%g_hat(2, ih)) cycle
1825 IF (pw_grid%g_hat(3, ig) /= ref_grid%g_hat(3, ih)) cycle
1826 pw_grid%gidx(ig) = ih
1827 EXIT
1828 END DO
1829 IF (pw_grid%gidx(ig) == 0) THEN
1830 WRITE (*, "(A,2I10)") " G-Shell ", is + 1, ip + 1
1831 WRITE (*, "(A,I10,3I6,F20.10)") &
1832 " G-vector", ig, pw_grid%g_hat(1:3, ig), pw_grid%gsq(ig)
1833 DO ih = 1, SIZE(ref_grid%gsq)
1834 IF (pw_grid%g_hat(1, ig) /= ref_grid%g_hat(1, ih)) cycle
1835 IF (pw_grid%g_hat(2, ig) /= ref_grid%g_hat(2, ih)) cycle
1836 IF (pw_grid%g_hat(3, ig) /= ref_grid%g_hat(3, ih)) cycle
1837 WRITE (*, "(A,I10,3I6,F20.10)") &
1838 " G-vector", ih, ref_grid%g_hat(1:3, ih), ref_grid%gsq(ih)
1839 END DO
1840 cpabort("G vector not found")
1841 END IF
1842 is = ip
1843 END DO
1844 ELSE
1845 ! for the case pw_grid > ref_grid
1846 is = it
1847 DO ig = it + 1, ngr
1848 gig = ref_grid%gsq(ig)
1849 gigr = max(1._dp, gig)
1850 g_found = .false.
1851 DO ih = is + 1, ng
1852 IF (abs(pw_grid%gsq(ih) - gig)/gigr > 1.e-12_dp) cycle
1853 g_found = .true.
1854 EXIT
1855 END DO
1856 IF (.NOT. g_found) THEN
1857 WRITE (*, "(A,I10,F20.10)") "G-vector", ig, ref_grid%gsq(ig)
1858 cpabort("G vector not found")
1859 END IF
1860 ip = ih - 1
1861 DO ih = ip + 1, ng
1862 IF (abs(pw_grid%gsq(ih) - gig)/gigr > 1.e-12_dp) cycle
1863 IF (pw_grid%g_hat(1, ih) /= ref_grid%g_hat(1, ig)) cycle
1864 IF (pw_grid%g_hat(2, ih) /= ref_grid%g_hat(2, ig)) cycle
1865 IF (pw_grid%g_hat(3, ih) /= ref_grid%g_hat(3, ig)) cycle
1866 pw_grid%gidx(ig) = ih
1867 EXIT
1868 END DO
1869 IF (pw_grid%gidx(ig) == 0) THEN
1870 WRITE (*, "(A,2I10)") " G-Shell ", is + 1, ip + 1
1871 WRITE (*, "(A,I10,3I6,F20.10)") &
1872 " G-vector", ig, ref_grid%g_hat(1:3, ig), ref_grid%gsq(ig)
1873 DO ih = 1, ng
1874 IF (pw_grid%g_hat(1, ih) /= ref_grid%g_hat(1, ig)) cycle
1875 IF (pw_grid%g_hat(2, ih) /= ref_grid%g_hat(2, ig)) cycle
1876 IF (pw_grid%g_hat(3, ih) /= ref_grid%g_hat(3, ig)) cycle
1877 WRITE (*, "(A,I10,3I6,F20.10)") &
1878 " G-vector", ih, pw_grid%g_hat(1:3, ih), pw_grid%gsq(ih)
1879 END DO
1880 cpabort("G vector not found")
1881 END IF
1882 is = ip
1883 END DO
1884 END IF
1885 ! test if all g-vectors are associated
1886 IF (any(pw_grid%gidx == 0)) THEN
1887 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 IF (pw_grid%have_g0) THEN
1894 IF (pw_grid%gsq(1) /= 0.0_dp) THEN
1895 cpabort("G = 0 not in first position")
1896 END IF
1897 END IF
1898
1899 CALL timestop(handle)
1900
1901 END SUBROUTINE pw_grid_sort
1902
1903! **************************************************************************************************
1904!> \brief ...
1905!> \param gsq ...
1906!> \param g_hat ...
1907!> \param idx ...
1908! **************************************************************************************************
1909 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 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 ng = SIZE(gsq)
1928 s_begin = -1.0_dp
1929 s1 = 0
1930 s2 = 0
1931 ig = 1
1932 DO ig = 1, ng
1933 IF (abs(gsq(ig) - s_begin) < small) THEN
1934 s2 = ig
1935 ELSE
1936 CALL redist(g_hat, idx, s1, s2)
1937 s_begin = gsq(ig)
1938 s1 = ig
1939 s2 = ig
1940 END IF
1941 END DO
1942 CALL redist(g_hat, idx, s1, s2)
1943
1944 CALL timestop(handle)
1945
1946 END SUBROUTINE sort_shells
1947
1948! **************************************************************************************************
1949!> \brief ...
1950!> \param g_hat ...
1951!> \param idx ...
1952!> \param s1 ...
1953!> \param s2 ...
1954! **************************************************************************************************
1955 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: indl
1964 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: slen
1965
1966 IF (s2 <= s1) RETURN
1967 ns = s2 - s1 + 1
1968 ALLOCATE (indl(ns))
1969 ALLOCATE (slen(ns))
1970
1971 DO i = s1, s2
1972 ii = idx(i)
1973 n1 = g_hat(1, ii)
1974 n2 = g_hat(2, ii)
1975 n3 = g_hat(3, ii)
1976 slen(i - s1 + 1) = 1000.0_dp*real(n1, dp) + &
1977 REAL(n2, dp) + 0.001_dp*REAL(n3, dp)
1978 END DO
1979 CALL sort(slen, ns, indl)
1980 DO i = 1, ns
1981 ii = indl(i) + s1 - 1
1982 indl(i) = idx(ii)
1983 END DO
1984 idx(s1:s2) = indl(1:ns)
1985
1986 DEALLOCATE (indl)
1987 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 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 IF (pw_grid%para%mode == pw_mode_local) RETURN
2010
2011 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 yz = 0
2017 pw_grid%para%yzp = 0
2018 pw_grid%para%yzq = 0
2019
2020 DO gpt = 1, SIZE(pw_grid%gsq)
2021 m = posm(pw_grid%g_hat(2, gpt)) + 1
2022 n = posn(pw_grid%g_hat(3, gpt)) + 1
2023 yz(m, n) = yz(m, n) + 1
2024 END DO
2025 IF (pw_grid%grid_span == halfspace) THEN
2026 DO gpt = 1, SIZE(pw_grid%gsq)
2027 m = negm(pw_grid%g_hat(2, gpt)) + 1
2028 n = negn(pw_grid%g_hat(3, gpt)) + 1
2029 yz(m, n) = yz(m, n) + 1
2030 END DO
2031 END IF
2032
2033 ip = pw_grid%para%my_pos
2034 is = 0
2035 DO i = 1, nz
2036 DO j = 1, ny
2037 IF (yz(j, i) > 0) THEN
2038 is = is + 1
2039 pw_grid%para%yzp(1, is, ip) = j
2040 pw_grid%para%yzp(2, is, ip) = i
2041 pw_grid%para%yzq(j, i) = is
2042 END IF
2043 END DO
2044 END DO
2045 END associate
2046
2047 cpassert(is == pw_grid%para%nyzray(ip))
2048 CALL pw_grid%para%group%sum(pw_grid%para%yzp)
2049
2050 CALL timestop(handle)
2051
2052 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: g
2076
2077 cell_deth = abs(det_3x3(cell_hmat))
2078 IF (cell_deth < 1.0e-10_dp) THEN
2079 CALL cp_abort(__location__, &
2080 "An invalid set of cell vectors was specified. "// &
2081 "The determinant det(h) is too small")
2082 END IF
2083 cell_h_inv = inv_3x3(cell_hmat)
2084
2085 g => pw_grid%g
2086
2087 CALL cell2grid(cell_hmat, cell_h_inv, cell_deth, pw_grid)
2088
2089 DO gpt = 1, SIZE(g, 2)
2090
2091 l = twopi*real(pw_grid%g_hat(1, gpt), kind=dp)
2092 m = twopi*real(pw_grid%g_hat(2, gpt), kind=dp)
2093 n = twopi*real(pw_grid%g_hat(3, gpt), kind=dp)
2094
2095 g(1, gpt) = l*cell_h_inv(1, 1) + m*cell_h_inv(2, 1) + n*cell_h_inv(3, 1)
2096 g(2, gpt) = l*cell_h_inv(1, 2) + m*cell_h_inv(2, 2) + n*cell_h_inv(3, 2)
2097 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 + g(3, gpt)*g(3, gpt)
2102
2103 END DO
2104
2105 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 SUBROUTINE pw_grid_retain(pw_grid)
2117 TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid
2118
2119 cpassert(pw_grid%ref_count > 0)
2120 pw_grid%ref_count = pw_grid%ref_count + 1
2121 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 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 IF (ASSOCIATED(pw_grid)) THEN
2142 cpassert(pw_grid%ref_count > 0)
2143 pw_grid%ref_count = pw_grid%ref_count - 1
2144 IF (pw_grid%ref_count == 0) THEN
2145 IF (ASSOCIATED(pw_grid%gidx)) THEN
2146 DEALLOCATE (pw_grid%gidx)
2147 END IF
2148 IF (ASSOCIATED(pw_grid%g)) THEN
2149 DEALLOCATE (pw_grid%g)
2150 END IF
2151 IF (ASSOCIATED(pw_grid%gsq)) THEN
2152 DEALLOCATE (pw_grid%gsq)
2153 END IF
2154 IF (ALLOCATED(pw_grid%g_hat)) THEN
2155 DEALLOCATE (pw_grid%g_hat)
2156 END IF
2157 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 DEALLOCATE (pw_grid%g_hatmap)
2164#endif
2165 END IF
2166 IF (ASSOCIATED(pw_grid%grays)) THEN
2167 DEALLOCATE (pw_grid%grays)
2168 END IF
2169 IF (ALLOCATED(pw_grid%mapl%pos)) THEN
2170 DEALLOCATE (pw_grid%mapl%pos)
2171 END IF
2172 IF (ALLOCATED(pw_grid%mapm%pos)) THEN
2173 DEALLOCATE (pw_grid%mapm%pos)
2174 END IF
2175 IF (ALLOCATED(pw_grid%mapn%pos)) THEN
2176 DEALLOCATE (pw_grid%mapn%pos)
2177 END IF
2178 IF (ALLOCATED(pw_grid%mapl%neg)) THEN
2179 DEALLOCATE (pw_grid%mapl%neg)
2180 END IF
2181 IF (ALLOCATED(pw_grid%mapm%neg)) THEN
2182 DEALLOCATE (pw_grid%mapm%neg)
2183 END IF
2184 IF (ALLOCATED(pw_grid%mapn%neg)) THEN
2185 DEALLOCATE (pw_grid%mapn%neg)
2186 END IF
2187 IF (ALLOCATED(pw_grid%para%bo)) THEN
2188 DEALLOCATE (pw_grid%para%bo)
2189 END IF
2190 IF (pw_grid%para%mode == pw_mode_distributed) THEN
2191 IF (ALLOCATED(pw_grid%para%yzp)) THEN
2192 DEALLOCATE (pw_grid%para%yzp)
2193 END IF
2194 IF (ALLOCATED(pw_grid%para%yzq)) THEN
2195 DEALLOCATE (pw_grid%para%yzq)
2196 END IF
2197 IF (ALLOCATED(pw_grid%para%nyzray)) THEN
2198 DEALLOCATE (pw_grid%para%nyzray)
2199 END IF
2200 END IF
2201 ! also release groups
2202 CALL pw_grid%para%group%free()
2203 IF (product(pw_grid%para%rs_dims) /= 0) &
2204 CALL pw_grid%para%rs_group%free()
2205 IF (ALLOCATED(pw_grid%para%pos_of_x)) THEN
2206 DEALLOCATE (pw_grid%para%pos_of_x)
2207 END IF
2208
2209 IF (ASSOCIATED(pw_grid)) THEN
2210 DEALLOCATE (pw_grid)
2211 END IF
2212 END IF
2213 END IF
2214 NULLIFY (pw_grid)
2215 END SUBROUTINE pw_grid_release
2216
2217END MODULE pw_grids
idx
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition grid_common.h:153
mathlib::det_3x3
Definition mathlib.F:65
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::int_size
integer, parameter, public int_size
Definition kinds.F:36
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::inv_3x3
pure real(kind=dp) function, dimension(3, 3), public inv_3x3(a)
Returns the inverse of the 3 x 3 matrix a.
Definition mathlib.F:516
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
message_passing::mp_dims_create
subroutine, public mp_dims_create(nodes, dims)
wrapper to MPI_Dims_create
Definition message_passing.F:2132
message_passing::mp_comm_self
type(mp_comm_type), parameter, public mp_comm_self
Definition message_passing.F:812
offload_api
Fortran API for the offload package, which is written in C.
Definition offload_api.F:12
offload_api::offload_activate_chosen_device
subroutine, public offload_activate_chosen_device()
Activates the device selected via offload_set_chosen_device()
Definition offload_api.F:174
offload_api::offload_free_pinned_mem
integer function, public offload_free_pinned_mem(buffer)
free pinned memory
Definition offload_api.F:75
offload_api::offload_malloc_pinned_mem
integer function, public offload_malloc_pinned_mem(buffer, length)
allocate pinned memory.
Definition offload_api.F:52
pw_grid_info
This module returns additional info on PW grids.
Definition pw_grid_info.F:14
pw_grid_info::pw_find_cutoff
real(kind=dp) function, public pw_find_cutoff(npts, h_inv)
Given a grid and a box, calculate the corresponding cutoff *** This routine calculates the cutoff in ...
Definition pw_grid_info.F:272
pw_grid_info::pw_grid_bounds_from_n
integer function, dimension(2, 3), public pw_grid_bounds_from_n(npts)
returns the bounds that distribute n points evenly around 0
Definition pw_grid_info.F:248
pw_grid_info::pw_grid_init_setup
integer function, dimension(3), public pw_grid_init_setup(hmat, cutoff, spherical, odd, fft_usage, ncommensurate, icommensurate, ref_grid, n_orig)
...
Definition pw_grid_info.F:48
pw_grid_types
Definition pw_grid_types.F:13
pw_grid_types::halfspace
integer, parameter, public halfspace
Definition pw_grid_types.F:28
pw_grid_types::pw_mode_local
integer, parameter, public pw_mode_local
Definition pw_grid_types.F:29
pw_grid_types::fullspace
integer, parameter, public fullspace
Definition pw_grid_types.F:28
pw_grid_types::pw_mode_distributed
integer, parameter, public pw_mode_distributed
Definition pw_grid_types.F:29
pw_grids
This module defines the grid data type and some basic operations on it.
Definition pw_grids.F:36
pw_grids::pw_grid_compare
logical function, public pw_grid_compare(grida, gridb)
Check if two pw_grids are equal.
Definition pw_grids.F:147
pw_grids::do_pw_grid_blocked_false
integer, parameter, public do_pw_grid_blocked_false
Definition pw_grids.F:78
pw_grids::pw_grid_release
subroutine, public pw_grid_release(pw_grid)
releases the given pw grid
Definition pw_grids.F:2133
pw_grids::do_pw_grid_blocked_true
integer, parameter, public do_pw_grid_blocked_true
Definition pw_grids.F:78
pw_grids::pw_grid_retain
subroutine, public pw_grid_retain(pw_grid)
retains the given pw grid
Definition pw_grids.F:2117
pw_grids::do_pw_grid_blocked_free
integer, parameter, public do_pw_grid_blocked_free
Definition pw_grids.F:78
pw_grids::pw_grid_setup
subroutine, public pw_grid_setup(cell_hmat, pw_grid, grid_span, cutoff, bounds, bounds_local, npts, spherical, odd, fft_usage, ncommensurate, icommensurate, blocked, ref_grid, rs_dims, iounit)
sets up a pw_grid
Definition pw_grids.F:286
pw_grids::pw_grid_create
subroutine, public pw_grid_create(pw_grid, pe_group, local)
Initialize a PW grid with all defaults.
Definition pw_grids.F:93
pw_grids::get_pw_grid_info
subroutine, public get_pw_grid_info(pw_grid, id_nr, mode, vol, dvol, npts, ngpts, ngpts_cut, dr, cutoff, orthorhombic, gvectors, gsquare)
Access to information stored in the pw_grid_type.
Definition pw_grids.F:184
pw_grids::pw_grid_change
subroutine, public pw_grid_change(cell_hmat, pw_grid)
Recalculate the g-vectors after a change of the box.
Definition pw_grids.F:2068
util
All kind of helpful little routines.
Definition util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition util.F:333
offload_activate_chosen_device
void offload_activate_chosen_device(void)
Activates the device selected via offload_set_chosen_device()
Definition offload_library.c:99
message_passing::mp_comm_type
Definition message_passing.F:189
pw_grid_types::map_pn
Definition pw_grid_types.F:33
pw_grid_types::pw_grid_type
Definition pw_grid_types.F:62