LCOV - code coverage report
Current view: top level - src/pw - pw_grids.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:58e3e09) Lines: 857 932 92.0 %
Date: 2024-03-29 07:50:05 Functions: 23 23 100.0 %

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

Generated by: LCOV version 1.15