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