(git:b279b6b)
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 ! **************************************************************************************************
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  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 
2217 END MODULE pw_grids
det_3x3
real(kind=dp) function det_3x3(a)
...
Definition: dumpdcd.F:1091
idx
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition: grid_common.h:153
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