LCOV - code coverage report
Current view: top level - src/pw - pw_grids.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 91.9 % 938 862
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 23 23

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

Generated by: LCOV version 2.0-1