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