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