|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00034 MODULE pw_grids 00035 USE cell_types, ONLY: cell_type 00036 USE f77_blas 00037 USE kinds, ONLY: dp,& 00038 int_8 00039 USE mathconstants, ONLY: twopi 00040 USE message_passing, ONLY: & 00041 MPI_COMM_SELF, mp_cart_coords, mp_cart_create, mp_cart_rank, & 00042 mp_comm_compare, mp_comm_dup, mp_comm_free, mp_dims_create, & 00043 mp_environ, mp_max, mp_min, mp_sum 00044 USE pw_grid_info, ONLY: pw_find_cutoff,& 00045 pw_grid_bounds_from_n,& 00046 pw_grid_init_setup 00047 USE pw_grid_types, ONLY: & 00048 FULLSPACE, HALFSPACE, PW_MODE_DISTRIBUTED, PW_MODE_LOCAL, & 00049 do_pw_grid_blocked_false, do_pw_grid_blocked_free, & 00050 do_pw_grid_blocked_true, map_pn, pw_grid_type 00051 USE termination, ONLY: stop_program 00052 USE timings, ONLY: timeset,& 00053 timestop 00054 USE util, ONLY: get_limit,& 00055 sort 00056 #include "cp_common_uses.h" 00057 00058 IMPLICIT NONE 00059 00060 PRIVATE 00061 PUBLIC :: pw_grid_create, pw_grid_retain, pw_grid_release 00062 PUBLIC :: get_pw_grid_info, set_pw_grid_info, pw_grid_compare 00063 PUBLIC :: pw_grid_setup 00064 PUBLIC :: pw_grid_change, create_gvectors 00065 00066 INTEGER :: grid_tag = 0 00067 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_grids' 00068 00069 CONTAINS 00070 00071 ! ***************************************************************************** 00077 SUBROUTINE pw_grid_create ( pw_grid, pe_group, local, error ) 00078 00079 TYPE(pw_grid_type), POINTER :: pw_grid 00080 INTEGER, INTENT(in) :: pe_group 00081 LOGICAL, INTENT(IN), OPTIONAL :: local 00082 TYPE(cp_error_type), INTENT(inout) :: error 00083 00084 CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_create', 00085 routineP = moduleN//':'//routineN 00086 00087 INTEGER :: stat 00088 LOGICAL :: failure, my_local 00089 00090 failure =.FALSE. 00091 my_local=.FALSE. 00092 IF (PRESENT(local)) my_local = local 00093 CPPrecondition(.NOT.ASSOCIATED(pw_grid),cp_failure_level,routineP,error,failure) 00094 IF (.NOT.failure) THEN 00095 ALLOCATE(pw_grid,stat=stat) 00096 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 00097 END IF 00098 IF (.NOT.failure) THEN 00099 pw_grid % bounds = 0 00100 pw_grid % cutoff = 0.0_dp 00101 pw_grid % grid_span = FULLSPACE 00102 pw_grid % para % mode = PW_MODE_LOCAL 00103 pw_grid % para % rs_dims = 0 00104 pw_grid % reference = 0 00105 pw_grid % ref_count = 1 00106 NULLIFY ( pw_grid % g ) 00107 NULLIFY ( pw_grid % gsq ) 00108 NULLIFY ( pw_grid % g_hat ) 00109 NULLIFY ( pw_grid % gidx ) 00110 NULLIFY ( pw_grid % grays ) 00111 NULLIFY ( pw_grid % mapl % pos ) 00112 NULLIFY ( pw_grid % mapl % neg ) 00113 NULLIFY ( pw_grid % mapm % pos ) 00114 NULLIFY ( pw_grid % mapm % neg ) 00115 NULLIFY ( pw_grid % mapn % pos ) 00116 NULLIFY ( pw_grid % mapn % neg ) 00117 NULLIFY ( pw_grid % para % yzp ) 00118 NULLIFY ( pw_grid % para % yzq ) 00119 NULLIFY ( pw_grid % para % nyzray ) 00120 NULLIFY ( pw_grid % para % bo ) 00121 NULLIFY ( pw_grid % para % pos_of_x ) 00122 00123 ! assign a unique tag to this grid 00124 grid_tag = grid_tag + 1 00125 pw_grid %id_nr = grid_tag 00126 00127 ! parallel info 00128 CALL mp_comm_dup ( pe_group, pw_grid % para % group ) 00129 CALL mp_environ ( pw_grid % para % group_size, & 00130 pw_grid % para % my_pos, & 00131 pw_grid % para % group ) 00132 pw_grid % para % group_head_id = 0 00133 pw_grid % para % group_head = & 00134 ( pw_grid % para % group_head_id == pw_grid % para % my_pos ) 00135 IF (pw_grid % para % group_size > 1 .AND. (.NOT.my_local)) THEN 00136 pw_grid % para % mode = PW_MODE_DISTRIBUTED 00137 ELSE 00138 pw_grid % para % mode = PW_MODE_LOCAL 00139 END IF 00140 END IF 00141 00142 END SUBROUTINE pw_grid_create 00143 00144 ! ***************************************************************************** 00150 FUNCTION pw_grid_compare ( grida, gridb ) RESULT ( equal ) 00151 00152 TYPE(pw_grid_type), INTENT(IN) :: grida, gridb 00153 LOGICAL :: equal 00154 00155 !------------------------------------------------------------------------------ 00156 00157 IF ( grida %id_nr == gridb %id_nr ) THEN 00158 equal = .TRUE. 00159 ELSE 00160 ! for the moment all grids with different identifiers are considered as not equal 00161 ! later we can get this more relaxed 00162 equal = .FALSE. 00163 END IF 00164 00165 END FUNCTION pw_grid_compare 00166 00167 ! ***************************************************************************** 00173 SUBROUTINE get_pw_grid_info(pw_grid, id_nr, mode, vol, dvol, npts, ngpts,& 00174 ngpts_cut,dr, cutoff, orthorhombic, gvectors,gsquare,error) 00175 00176 TYPE(pw_grid_type), POINTER :: pw_grid 00177 INTEGER, INTENT(OUT), OPTIONAL :: id_nr, mode 00178 REAL(dp), INTENT(OUT), OPTIONAL :: vol, dvol 00179 INTEGER, DIMENSION(3), INTENT(OUT), 00180 OPTIONAL :: npts 00181 INTEGER(int_8), INTENT(OUT), OPTIONAL :: ngpts, ngpts_cut 00182 REAL(dp), DIMENSION(3), INTENT(OUT), 00183 OPTIONAL :: dr 00184 REAL(dp), INTENT(OUT), OPTIONAL :: cutoff 00185 LOGICAL, INTENT(OUT), OPTIONAL :: orthorhombic 00186 REAL(dp), DIMENSION(:, :), OPTIONAL, 00187 POINTER :: gvectors 00188 REAL(dp), DIMENSION(:), OPTIONAL, 00189 POINTER :: gsquare 00190 TYPE(cp_error_type), INTENT(INOUT) :: error 00191 00192 CHARACTER(len=*), PARAMETER :: routineN = 'get_pw_grid_info', 00193 routineP = moduleN//':'//routineN 00194 00195 LOGICAL :: failure 00196 00197 failure=.FALSE. 00198 CPPrecondition(pw_grid%ref_count>0,cp_failure_level,routineP,error,failure) 00199 00200 IF ( .NOT. failure ) THEN 00201 IF ( PRESENT(id_nr) ) id_nr = pw_grid%id_nr 00202 IF ( PRESENT(mode) ) mode = pw_grid%para%mode 00203 IF ( PRESENT(vol) ) vol = pw_grid%vol 00204 IF ( PRESENT(dvol) ) dvol = pw_grid%dvol 00205 IF ( PRESENT(npts) ) npts(1:3) = pw_grid%npts(1:3) 00206 IF ( PRESENT(ngpts) ) ngpts = pw_grid%ngpts 00207 IF ( PRESENT(ngpts_cut) ) ngpts_cut = pw_grid%ngpts_cut 00208 IF ( PRESENT(dr) ) dr = pw_grid%dr 00209 IF ( PRESENT(cutoff) ) cutoff = pw_grid%cutoff 00210 IF ( PRESENT(orthorhombic) ) orthorhombic = pw_grid%orthorhombic 00211 IF ( PRESENT(gvectors) ) gvectors => pw_grid%g 00212 IF ( PRESENT(gsquare) ) gsquare => pw_grid%gsq 00213 END IF 00214 00215 END SUBROUTINE get_pw_grid_info 00216 00217 ! ***************************************************************************** 00223 SUBROUTINE set_pw_grid_info(pw_grid, grid_span, npts, bounds, cutoff, spherical,& 00224 error) 00225 00226 TYPE(pw_grid_type), POINTER :: pw_grid 00227 INTEGER, INTENT(in), OPTIONAL :: grid_span 00228 INTEGER, DIMENSION(3), INTENT(IN), 00229 OPTIONAL :: npts 00230 INTEGER, DIMENSION(2, 3), INTENT(IN), 00231 OPTIONAL :: bounds 00232 REAL(KIND=dp), INTENT(IN), OPTIONAL :: cutoff 00233 LOGICAL, INTENT(IN), OPTIONAL :: spherical 00234 TYPE(cp_error_type), INTENT(INOUT) :: error 00235 00236 CHARACTER(len=*), PARAMETER :: routineN = 'set_pw_grid_info', 00237 routineP = moduleN//':'//routineN 00238 00239 LOGICAL :: failure 00240 00241 failure=.FALSE. 00242 CPPrecondition(pw_grid%ref_count>0,cp_failure_level,routineP,error,failure) 00243 00244 IF ( .NOT. failure ) THEN 00245 IF ( PRESENT(grid_span) ) THEN 00246 pw_grid%grid_span = grid_span 00247 END IF 00248 IF ( PRESENT(bounds) .AND. PRESENT(npts) ) THEN 00249 pw_grid%bounds = bounds 00250 pw_grid%npts = npts 00251 CPPostcondition(ALL(npts==bounds(2,:)-bounds(1,:)+1),cp_failure_level,routineP,error,failure) 00252 ELSE IF ( PRESENT(bounds) ) THEN 00253 pw_grid%bounds = bounds 00254 pw_grid%npts = bounds(2,:)-bounds(1,:)+1 00255 ELSE IF ( PRESENT(npts) ) THEN 00256 pw_grid%npts = npts 00257 pw_grid%bounds = pw_grid_bounds_from_n(npts,error) 00258 END IF 00259 IF ( PRESENT(cutoff) ) THEN 00260 pw_grid%cutoff = cutoff 00261 IF( PRESENT(spherical) ) THEN 00262 pw_grid%spherical = spherical 00263 ELSE 00264 pw_grid%spherical = .FALSE. 00265 END IF 00266 END IF 00267 END IF 00268 00269 END SUBROUTINE set_pw_grid_info 00270 00271 ! ***************************************************************************** 00279 SUBROUTINE pw_grid_setup ( cell, pw_grid, grid_span, cutoff, bounds, npts, & 00280 spherical, odd, fft_usage, ncommensurate, icommensurate, blocked, ref_grid,& 00281 rs_dims, iounit, error ) 00282 00283 TYPE(cell_type), INTENT(IN) :: cell 00284 TYPE(pw_grid_type), POINTER :: pw_grid 00285 INTEGER, INTENT(in), OPTIONAL :: grid_span 00286 REAL(KIND=dp), INTENT(IN), OPTIONAL :: cutoff 00287 INTEGER, DIMENSION(2, 3), INTENT(IN), 00288 OPTIONAL :: bounds 00289 INTEGER, DIMENSION(3), INTENT(IN), 00290 OPTIONAL :: npts 00291 LOGICAL, INTENT(in), OPTIONAL :: spherical, odd, fft_usage 00292 INTEGER, INTENT(in), OPTIONAL :: ncommensurate, icommensurate, 00293 blocked 00294 TYPE(pw_grid_type), INTENT(IN), OPTIONAL :: ref_grid 00295 INTEGER, DIMENSION(2), INTENT(in), 00296 OPTIONAL :: rs_dims 00297 INTEGER, INTENT(in), OPTIONAL :: iounit 00298 TYPE(cp_error_type), INTENT(inout) :: error 00299 00300 CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_setup', 00301 routineP = moduleN//':'//routineN 00302 00303 INTEGER :: handle, my_icommensurate, 00304 my_ncommensurate 00305 INTEGER, DIMENSION(3) :: n 00306 LOGICAL :: my_fft_usage, my_odd, 00307 my_spherical 00308 REAL(KIND=dp) :: my_cutoff 00309 00310 CALL timeset(routineN,handle) 00311 00312 IF ( PRESENT(grid_span) ) THEN 00313 CALL set_pw_grid_info ( pw_grid, grid_span=grid_span, error=error) 00314 END IF 00315 00316 IF ( PRESENT(spherical) ) THEN 00317 my_spherical = spherical 00318 ELSE 00319 my_spherical = .FALSE. 00320 END IF 00321 00322 IF ( PRESENT(odd) ) THEN 00323 my_odd = odd 00324 ELSE 00325 my_odd = .FALSE. 00326 END IF 00327 00328 IF ( PRESENT(fft_usage) ) THEN 00329 my_fft_usage = fft_usage 00330 ELSE 00331 my_fft_usage = .FALSE. 00332 END IF 00333 00334 IF ( PRESENT(ncommensurate) ) THEN 00335 my_ncommensurate = ncommensurate 00336 IF ( PRESENT(icommensurate) ) THEN 00337 my_icommensurate = icommensurate 00338 ELSE 00339 my_icommensurate = MIN(1,ncommensurate) 00340 END IF 00341 ELSE 00342 my_ncommensurate = 0 00343 my_icommensurate = 1 00344 END IF 00345 00346 IF ( PRESENT(bounds) ) THEN 00347 IF ( PRESENT(cutoff) ) THEN 00348 CALL set_pw_grid_info ( pw_grid, bounds=bounds, cutoff=cutoff, & 00349 spherical=my_spherical, error=error) 00350 ELSE 00351 n = bounds(2,:) - bounds(1,:) + 1 00352 my_cutoff = pw_find_cutoff ( n, cell%h_inv, error=error) 00353 my_cutoff = 0.5_dp * my_cutoff * my_cutoff 00354 CALL set_pw_grid_info ( pw_grid, bounds=bounds, cutoff=my_cutoff, & 00355 spherical=my_spherical, error=error) 00356 END IF 00357 ELSE IF ( PRESENT(npts) ) THEN 00358 n = npts 00359 IF ( PRESENT(cutoff) ) THEN 00360 my_cutoff=cutoff 00361 ELSE 00362 my_cutoff = pw_find_cutoff ( npts, cell%h_inv, error=error) 00363 my_cutoff = 0.5_dp * my_cutoff * my_cutoff 00364 END IF 00365 IF ( my_fft_usage ) THEN 00366 n = pw_grid_init_setup(cell%hmat, cutoff=my_cutoff,& 00367 spherical=my_spherical, odd=my_odd, fft_usage=my_fft_usage,& 00368 ncommensurate=my_ncommensurate, icommensurate=my_icommensurate,& 00369 ref_grid=ref_grid, n_orig=n, error=error) 00370 END IF 00371 CALL set_pw_grid_info ( pw_grid, npts=n, cutoff=my_cutoff, & 00372 spherical=my_spherical, error=error) 00373 ELSE IF ( PRESENT(cutoff) ) THEN 00374 n = pw_grid_init_setup(cell%hmat,cutoff=cutoff,& 00375 spherical=my_spherical,odd=my_odd,fft_usage=my_fft_usage,& 00376 ncommensurate=my_ncommensurate,icommensurate=my_icommensurate,& 00377 ref_grid=ref_grid,error=error) 00378 CALL set_pw_grid_info ( pw_grid, npts=n, cutoff=cutoff, & 00379 spherical=my_spherical, error=error) 00380 ELSE 00381 CALL stop_program(routineN,moduleN,__LINE__,& 00382 "BOUNDS, NPTS or CUTOFF have to be specified") 00383 END IF 00384 00385 CALL pw_grid_setup_internal ( cell, pw_grid, blocked, ref_grid, rs_dims, iounit, error ) 00386 CALL timestop(handle) 00387 00388 END SUBROUTINE pw_grid_setup 00389 00390 ! ***************************************************************************** 00409 SUBROUTINE pw_grid_setup_internal ( cell, pw_grid, blocked, ref_grid, rs_dims, iounit, error ) 00410 00411 TYPE(cell_type), INTENT(IN) :: cell 00412 TYPE(pw_grid_type), POINTER :: pw_grid 00413 INTEGER, INTENT(in), OPTIONAL :: blocked 00414 TYPE(pw_grid_type), INTENT(in), OPTIONAL :: ref_grid 00415 INTEGER, DIMENSION(2), INTENT(in), 00416 OPTIONAL :: rs_dims 00417 INTEGER, INTENT(in), OPTIONAL :: iounit 00418 TYPE(cp_error_type), INTENT(inout) :: error 00419 00420 CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_setup_internal', 00421 routineP = moduleN//':'//routineN 00422 00423 INTEGER :: allocstat, ires, n(3) 00424 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: yz_mask 00425 LOGICAL :: failure 00426 REAL(KIND=dp) :: ecut 00427 00428 !------------------------------------------------------------------------------ 00429 00430 failure=.FALSE. 00431 CPPrecondition(ASSOCIATED(pw_grid),cp_failure_level,routineP,error,failure) 00432 CPPrecondition(pw_grid%ref_count>0,cp_failure_level,routineP,error,failure) 00433 00434 00435 ! set pointer to possible reference grid 00436 IF ( PRESENT ( ref_grid ) ) THEN 00437 pw_grid % reference = ref_grid %id_nr 00438 END IF 00439 00440 IF (pw_grid%spherical) THEN 00441 ecut = pw_grid%cutoff 00442 ELSE 00443 ecut = 1.e10_dp 00444 END IF 00445 00446 n ( : ) = pw_grid % npts ( : ) 00447 00448 ! Find the number of grid points 00449 ! yz_mask counts the number of g-vectors orthogonal to the yz plane 00450 ! the indices in yz_mask are from -n/2 .. n/2 shifted by n/2 + 1 00451 ! these are not mapped indices ! 00452 ALLOCATE ( yz_mask ( n(2), n(3) ), STAT = allocstat ) 00453 CPPrecondition(allocstat==0,cp_failure_level,routineP,error,failure) 00454 CALL pw_grid_count ( cell % h_inv, pw_grid, ecut, yz_mask ) 00455 00456 ! Check if reference grid is compatible 00457 IF ( PRESENT ( ref_grid ) ) THEN 00458 CPPrecondition(pw_grid % para % mode == ref_grid % para % mode,cp_failure_level,routineP,error,failure) 00459 IF ( pw_grid % para % mode == PW_MODE_DISTRIBUTED ) THEN 00460 CALL mp_comm_compare(pw_grid % para % group, ref_grid % para % group, ires ) 00461 CPPrecondition(ires <= 3,cp_failure_level,routineP,error,failure) !FM make it >3 ? 00462 END IF 00463 CPPrecondition(pw_grid % grid_span == ref_grid % grid_span,cp_failure_level,routineP,error,failure) 00464 CPPrecondition(pw_grid % spherical .EQV. ref_grid % spherical,cp_failure_level,routineP,error,failure) 00465 END IF 00466 00467 ! Distribute grid 00468 CALL pw_grid_distribute ( pw_grid, yz_mask, ref_grid, blocked, & 00469 rs_dims=rs_dims ,error=error) 00470 00471 ! Allocate the grid fields 00472 CALL pw_grid_allocate ( pw_grid, pw_grid % ngpts_cut_local, & 00473 pw_grid % bounds, error ) 00474 00475 ! Fill in the grid structure 00476 CALL pw_grid_assign ( cell % h_inv, pw_grid, ecut, error ) 00477 00478 ! Sort g vector wrt length (only local for each processor) 00479 CALL pw_grid_sort ( pw_grid, ref_grid, error ) 00480 00481 CALL pw_grid_remap ( pw_grid, yz_mask, error ) 00482 00483 DEALLOCATE ( yz_mask , STAT=allocstat ) 00484 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 00485 00486 pw_grid % vol = ABS ( cell % deth ) 00487 pw_grid % dvol = pw_grid % vol / REAL ( pw_grid % ngpts,KIND=dp) 00488 pw_grid % dr ( 1 ) = SQRT ( SUM ( cell % hmat ( :, 1 ) ** 2 ) ) & 00489 / REAL ( pw_grid % npts ( 1 ),KIND=dp) 00490 pw_grid % dr ( 2 ) = SQRT ( SUM ( cell % hmat ( :, 2 ) ** 2 ) ) & 00491 / REAL ( pw_grid % npts ( 2 ),KIND=dp) 00492 pw_grid % dr ( 3 ) = SQRT ( SUM ( cell % hmat ( :, 3 ) ** 2 ) ) & 00493 / REAL ( pw_grid % npts ( 3 ),KIND=dp) 00494 pw_grid % dh ( :,1 ) = cell % hmat ( :, 1 ) / REAL ( pw_grid % npts ( 1 ),KIND=dp) 00495 pw_grid % dh ( :,2 ) = cell % hmat ( :, 2 ) / REAL ( pw_grid % npts ( 2 ),KIND=dp) 00496 pw_grid % dh ( :,3 ) = cell % hmat ( :, 3 ) / REAL ( pw_grid % npts ( 3 ),KIND=dp) 00497 pw_grid % dh_inv ( 1,: ) = cell % h_inv ( 1,: ) * REAL ( pw_grid % npts ( 1 ),KIND=dp) 00498 pw_grid % dh_inv ( 2,: ) = cell % h_inv ( 2,: ) * REAL ( pw_grid % npts ( 2 ),KIND=dp) 00499 pw_grid % dh_inv ( 3,: ) = cell % h_inv ( 3,: ) * REAL ( pw_grid % npts ( 3 ),KIND=dp) 00500 pw_grid % orthorhombic = cell % orthorhombic 00501 00502 ! 00503 ! Output: All the information of this grid type 00504 ! 00505 00506 IF(PRESENT(iounit)) THEN 00507 CALL pw_grid_print ( pw_grid, iounit, error ) 00508 END IF 00509 00510 END SUBROUTINE pw_grid_setup_internal 00511 00512 ! ***************************************************************************** 00513 SUBROUTINE create_gvectors(pw_grid,cell,ecut,blocked,ref_grid,error) 00514 ! Find the number of grid points 00515 ! yz_mask counts the number of g-vectors orthogonal to the yz plane 00516 ! the indices in yz_mask are from -n/2 .. n/2 shifted by n/2 + 1 00517 ! these are not mapped indices ! 00518 00519 TYPE(pw_grid_type), POINTER :: pw_grid 00520 TYPE(cell_type), INTENT(IN) :: cell 00521 REAL(KIND=dp) :: ecut 00522 INTEGER, OPTIONAL :: blocked 00523 TYPE(pw_grid_type), INTENT(IN), OPTIONAL :: ref_grid 00524 TYPE(cp_error_type), INTENT(inout) :: error 00525 00526 CHARACTER(len=*), PARAMETER :: routineN = 'create_gvectors', 00527 routineP = moduleN//':'//routineN 00528 00529 INTEGER :: istat, n(3) 00530 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: yz_mask 00531 LOGICAL :: failure 00532 00533 failure = .FALSE. 00534 00535 n ( : ) = pw_grid % npts ( : ) 00536 00537 ALLOCATE ( yz_mask ( n(2), n(3) ), STAT = istat ) 00538 CPPrecondition(istat == 0,cp_failure_level,routineP,error,failure) 00539 CALL pw_grid_count ( cell % h_inv, pw_grid, ecut, yz_mask ) 00540 00541 ! Distribute grid 00542 CALL pw_grid_distribute ( pw_grid, yz_mask, ref_grid, blocked ,error=error) 00543 00544 ! Allocate the grid fields 00545 CALL pw_grid_allocate ( pw_grid, pw_grid % ngpts_cut_local, & 00546 pw_grid % bounds, error ) 00547 00548 ! Fill in the grid structure 00549 CALL pw_grid_assign ( cell % h_inv, pw_grid, ecut, error ) 00550 00551 ! Sort g vector wrt length (only local for each processor) 00552 CALL pw_grid_sort ( pw_grid, ref_grid, error ) 00553 00554 CALL pw_grid_remap ( pw_grid, yz_mask, error ) 00555 00556 DEALLOCATE ( yz_mask, STAT=istat ) 00557 CPPrecondition(istat == 0,cp_failure_level,routineP,error,failure) 00558 00559 END SUBROUTINE create_gvectors 00560 00561 ! ***************************************************************************** 00565 SUBROUTINE pw_grid_print ( pw_grid, info, error ) 00566 00567 TYPE(pw_grid_type), POINTER :: pw_grid 00568 INTEGER, INTENT(IN) :: info 00569 TYPE(cp_error_type), INTENT(inout) :: error 00570 00571 CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_print', 00572 routineP = moduleN//':'//routineN 00573 00574 INTEGER :: i, n(3) 00575 REAL(KIND=dp) :: rv(3,3) 00576 00577 !------------------------------------------------------------------------------ 00578 ! 00579 ! Output: All the information of this grid type 00580 ! 00581 00582 IF ( pw_grid % para % mode == PW_MODE_LOCAL) THEN 00583 IF (info >= 0 ) THEN 00584 WRITE ( info, '(/,A,T71,I10)' ) & 00585 " PW_GRID: Information for grid number ", pw_grid %id_nr 00586 WRITE ( info, '(" PW_GRID: Cutoff [a.u.]",T71,f10.1)' ) pw_grid % cutoff 00587 IF ( pw_grid % spherical ) THEN 00588 WRITE ( info, '(A,T78,A)' ) " PW_GRID: spherical cutoff: ","YES" 00589 WRITE ( info, '(A,T71,I10)' ) " PW_GRID: Grid points within cutoff", & 00590 pw_grid % ngpts_cut 00591 ELSE 00592 WRITE ( info, '(A,T78,A)' ) " PW_GRID: spherical cutoff: "," NO" 00593 END IF 00594 DO i = 1, 3 00595 WRITE ( info, '(A,I3,T30,2I8,T62,A,T71,I10)' ) " PW_GRID: Bounds ", & 00596 i, pw_grid % bounds ( 1, I ), pw_grid % bounds ( 2, I ), & 00597 "Points:",pw_grid % npts ( I ) 00598 END DO 00599 WRITE ( info, '(A,G12.4,T50,A,T67,F14.4)' ) & 00600 " PW_GRID: Volume element (a.u.^3)", & 00601 pw_grid % dvol," Volume (a.u.^3) :",pw_grid % vol 00602 IF ( pw_grid % grid_span == HALFSPACE ) THEN 00603 WRITE ( info, '(A,T72,A)' ) " PW_GRID: Grid span","HALFSPACE" 00604 ELSE 00605 WRITE ( info, '(A,T72,A)' ) " PW_GRID: Grid span","FULLSPACE" 00606 END IF 00607 END IF 00608 ELSE 00609 00610 n ( 1 ) = pw_grid % ngpts_cut_local 00611 n ( 2 ) = pw_grid % ngpts_local 00612 CALL mp_sum ( n(1:2), pw_grid % para % group ) 00613 n ( 3 ) = SUM ( pw_grid % para % nyzray ) 00614 rv ( :, 1 ) = REAL ( n,KIND=dp) / REAL ( pw_grid % para % group_size,KIND=dp) 00615 n ( 1 ) = pw_grid % ngpts_cut_local 00616 n ( 2 ) = pw_grid % ngpts_local 00617 CALL mp_max ( n(1:2), pw_grid % para % group ) 00618 n ( 3 ) = MAXVAL ( pw_grid % para % nyzray ) 00619 rv ( :, 2 ) = REAL ( n,KIND=dp) 00620 n ( 1 ) = pw_grid % ngpts_cut_local 00621 n ( 2 ) = pw_grid % ngpts_local 00622 CALL mp_min ( n(1:2), pw_grid % para % group ) 00623 n ( 3 ) = MINVAL ( pw_grid % para % nyzray ) 00624 rv ( :, 3 ) = REAL ( n,KIND=dp) 00625 00626 IF ( pw_grid % para % group_head .AND. info>=0) THEN 00627 WRITE ( info, '(/,A,T71,I10)' ) & 00628 " PW_GRID: Information for grid number ", pw_grid %id_nr 00629 WRITE ( info, '(A,T60,I10,A)' ) & 00630 " PW_GRID: Grid distributed over ", pw_grid % para % group_size, & 00631 " processors" 00632 WRITE ( info, '(A,T71,2I5)' ) & 00633 " PW_GRID: Real space group dimensions ", pw_grid % para % rs_dims 00634 IF ( pw_grid % para % blocked ) THEN 00635 WRITE ( info, '(A,T78,A)' ) " PW_GRID: the grid is blocked: ","YES" 00636 ELSE 00637 WRITE ( info, '(A,T78,A)' ) " PW_GRID: the grid is blocked: "," NO" 00638 END IF 00639 WRITE ( info, '(" PW_GRID: Cutoff [a.u.]",T71,f10.1)' ) pw_grid % cutoff 00640 IF ( pw_grid % spherical ) THEN 00641 WRITE ( info, '(A,T78,A)' ) " PW_GRID: spherical cutoff: ","YES" 00642 WRITE ( info, '(A,T71,I10)' ) " PW_GRID: Grid points within cutoff", & 00643 pw_grid % ngpts_cut 00644 ELSE 00645 WRITE ( info, '(A,T78,A)' ) " PW_GRID: spherical cutoff: "," NO" 00646 END IF 00647 DO i = 1, 3 00648 WRITE ( info, '(A,I3,T30,2I8,T62,A,T71,I10)' ) " PW_GRID: Bounds ", & 00649 i, pw_grid % bounds ( 1, I ), pw_grid % bounds ( 2, I ), & 00650 "Points:",pw_grid % npts ( I ) 00651 END DO 00652 WRITE ( info, '(A,G12.4,T50,A,T67,F14.4)' ) & 00653 " PW_GRID: Volume element (a.u.^3)", & 00654 pw_grid % dvol," Volume (a.u.^3) :",pw_grid % vol 00655 IF ( pw_grid % grid_span == HALFSPACE ) THEN 00656 WRITE ( info, '(A,T72,A)' ) " PW_GRID: Grid span","HALFSPACE" 00657 ELSE 00658 WRITE ( info, '(A,T72,A)' ) " PW_GRID: Grid span","FULLSPACE" 00659 END IF 00660 WRITE ( info, '(A,T48,A)' ) " PW_GRID: Distribution", & 00661 " Average Max Min" 00662 WRITE ( info, '(A,T45,F12.1,2I12)' ) " PW_GRID: G-Vectors", & 00663 rv ( 1, 1 ), NINT ( rv ( 1, 2 ) ), NINT ( rv ( 1, 3 ) ) 00664 WRITE ( info, '(A,T45,F12.1,2I12)' ) " PW_GRID: G-Rays ", & 00665 rv ( 3, 1 ), NINT ( rv ( 3, 2 ) ), NINT ( rv ( 3, 3 ) ) 00666 WRITE ( info, '(A,T45,F12.1,2I12)' ) " PW_GRID: Real Space Points", & 00667 rv ( 2, 1 ), NINT ( rv ( 2, 2 ) ), NINT ( rv ( 2, 3 ) ) 00668 END IF ! group head 00669 END IF ! local 00670 00671 END SUBROUTINE pw_grid_print 00672 00673 ! ***************************************************************************** 00681 SUBROUTINE pw_grid_distribute ( pw_grid, yz_mask, ref_grid, blocked, rs_dims, error) 00682 00683 TYPE(pw_grid_type), POINTER :: pw_grid 00684 INTEGER, DIMENSION(:, :), INTENT(INOUT) :: yz_mask 00685 TYPE(pw_grid_type), INTENT(IN), OPTIONAL :: ref_grid 00686 INTEGER, INTENT(IN), OPTIONAL :: blocked 00687 INTEGER, DIMENSION(2), INTENT(in), 00688 OPTIONAL :: rs_dims 00689 TYPE(cp_error_type), INTENT(inout) :: error 00690 00691 CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_distribute', 00692 routineP = moduleN//':'//routineN 00693 00694 INTEGER :: blocked_local, coor( 2 ), gmax, handle, i, i1, i2, ierr, ip, 00695 ipl, ipp, itmp, j, l, lby, lbz, lo( 2 ), m, n, np, ns, nx, ny, nz, 00696 rsd( 2 ) 00697 INTEGER, ALLOCATABLE, DIMENSION(:) :: pemap 00698 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: yz_index 00699 LOGICAL :: blocking, failure 00700 00701 !------------------------------------------------------------------------------ 00702 00703 CALL timeset(routineN,handle) 00704 00705 lby = pw_grid % bounds ( 1, 2 ) 00706 lbz = pw_grid % bounds ( 1, 3 ) 00707 00708 pw_grid % ngpts = PRODUCT ( INT(pw_grid % npts,KIND=int_8) ) 00709 CPPrecondition(ALL(pw_grid%para%rs_dims==0),cp_failure_level,routineP,error,failure) 00710 IF (PRESENT(rs_dims)) THEN 00711 pw_grid%para%rs_dims=rs_dims 00712 END IF 00713 00714 00715 IF (PRESENT(blocked)) THEN 00716 blocked_local=blocked 00717 ELSE 00718 blocked_local=do_pw_grid_blocked_free 00719 ENDIF 00720 00721 00722 pw_grid % para % blocked = .FALSE. 00723 00724 IF ( pw_grid % para % mode == PW_MODE_LOCAL) THEN 00725 00726 pw_grid % bounds_local = pw_grid % bounds 00727 pw_grid % npts_local = pw_grid % npts 00728 pw_grid % ngpts_cut_local = pw_grid % ngpts_cut 00729 pw_grid % ngpts_local = PRODUCT ( pw_grid % npts_local ) 00730 pw_grid % para % rs_dims=1 00731 CALL mp_cart_create ( MPI_COMM_SELF, 2, & 00732 pw_grid % para % rs_dims, & 00733 pw_grid % para % rs_pos, & 00734 pw_grid % para % rs_group ) 00735 CALL mp_cart_rank ( pw_grid % para % rs_group, & 00736 pw_grid % para % rs_pos, & 00737 pw_grid % para % rs_mpo ) 00738 00739 ELSE 00740 00741 !..find the real space distribution 00742 nx = pw_grid % npts ( 1 ) 00743 ny = pw_grid % npts ( 2 ) 00744 nz = pw_grid % npts ( 3 ) 00745 np = pw_grid % para % group_size 00746 00747 00748 ! The user can specify 2 strictly positive indices => specific layout 00749 ! 1 strictly positive index => the other is fixed by the number of CPUs 00750 ! 0 strictly positive indices => fully free distribution 00751 ! if fully free or the user request can not be fulfilled, we optimize heuristically ourselves by: 00752 ! 1) nx>np -> taking a plane distribution (/np,1/) 00753 ! 2) nx<np -> taking the most square distribution 00754 ! if blocking is free: 00755 ! 1) blocked=.FALSE. for plane distributions 00756 ! 2) blocked=.TRUE. for non-plane distributions 00757 IF ( ANY(pw_grid % para % rs_dims<=0) ) THEN 00758 IF (ALL(pw_grid % para % rs_dims<=0) ) THEN 00759 pw_grid % para % rs_dims = (/0,0/) 00760 ELSE 00761 IF (pw_grid % para % rs_dims(1)>0) THEN 00762 pw_grid % para % rs_dims(2)=np/pw_grid % para % rs_dims(1) 00763 ELSE 00764 pw_grid % para % rs_dims(1)=np/pw_grid % para % rs_dims(2) 00765 ENDIF 00766 ENDIF 00767 ENDIF 00768 ! reset if the distribution can not be fulfilled 00769 IF ( PRODUCT(pw_grid % para % rs_dims) .NE. np ) pw_grid % para % rs_dims= (/0,0/) 00770 ! reset if the distribution can not be dealt with (/1,np/) 00771 IF ( ALL(pw_grid % para % rs_dims==(/1,np/)) ) pw_grid % para % rs_dims= (/0,0/) 00772 00773 ! if (/0,0/) now, we can optimize it ourselves 00774 IF ( ALL(pw_grid % para % rs_dims==(/0,0/)) ) THEN 00775 ! only small grids have a chance to be 2d distributed 00776 IF (nx<np) THEN 00777 ! gives the most square looking distribution 00778 CALL mp_dims_create ( np, pw_grid % para % rs_dims ) 00779 ! we tend to like the first index being smaller than the second 00780 IF (pw_grid % para % rs_dims(1)>pw_grid % para % rs_dims(2)) THEN 00781 itmp = pw_grid % para % rs_dims(1) 00782 pw_grid % para % rs_dims(1) = pw_grid % para % rs_dims(2) 00783 pw_grid % para % rs_dims(2) = itmp 00784 ENDIF 00785 ! but should avoid having the first index 1 in all cases 00786 IF (pw_grid % para % rs_dims(1) == 1) THEN 00787 itmp = pw_grid % para % rs_dims(1) 00788 pw_grid % para % rs_dims(1) = pw_grid % para % rs_dims(2) 00789 pw_grid % para % rs_dims(2) = itmp 00790 ENDIF 00791 ELSE 00792 pw_grid % para % rs_dims = (/np,1/) 00793 ENDIF 00794 ENDIF 00795 00796 ! now fix the blocking if we have a choice 00797 SELECT CASE(blocked_local) 00798 CASE(do_pw_grid_blocked_false) 00799 blocking=.FALSE. 00800 CASE(do_pw_grid_blocked_true) 00801 blocking=.TRUE. 00802 CASE(do_pw_grid_blocked_free) 00803 IF (ALL(pw_grid % para % rs_dims==(/np,1/))) THEN 00804 blocking=.FALSE. 00805 ELSE 00806 blocking=.TRUE. 00807 ENDIF 00808 CASE DEFAULT 00809 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00810 END SELECT 00811 00812 !..create group for real space distribution 00813 CALL mp_cart_create ( pw_grid % para % group, 2, & 00814 pw_grid % para % rs_dims, & 00815 pw_grid % para % rs_pos, & 00816 pw_grid % para % rs_group ) 00817 CALL mp_cart_rank ( pw_grid % para % rs_group, & 00818 pw_grid % para % rs_pos, & 00819 pw_grid % para % rs_mpo ) 00820 lo = get_limit ( nx, pw_grid % para % rs_dims ( 1 ), & 00821 pw_grid % para % rs_pos ( 1 ) ) 00822 pw_grid % bounds_local ( :, 1 ) = lo + pw_grid % bounds ( 1, 1 ) - 1 00823 lo = get_limit ( ny, pw_grid % para % rs_dims ( 2 ), & 00824 pw_grid % para % rs_pos ( 2 ) ) 00825 pw_grid % bounds_local ( :, 2 ) = lo + pw_grid % bounds ( 1, 2 ) - 1 00826 pw_grid % bounds_local ( :, 3 ) = pw_grid % bounds ( :, 3 ) 00827 pw_grid % npts_local ( : ) = pw_grid % bounds_local ( 2, : ) & 00828 - pw_grid % bounds_local ( 1, : ) + 1 00829 00830 !..the third distribution is needed for the second step in the FFT 00831 ALLOCATE ( pw_grid % para % bo ( 2, 3, 0:np-1, 3 ), STAT = ierr ) 00832 CPPrecondition(ierr == 0,cp_failure_level,routineP,error,failure) 00833 rsd = pw_grid % para % rs_dims 00834 DO ip = 0, np - 1 00835 CALL mp_cart_coords ( pw_grid % para % rs_group, ip, coor ) 00836 ! distribution xyZ 00837 pw_grid % para % bo ( 1:2, 1, ip, 1 ) = get_limit (nx,rsd(1),coor(1)) 00838 pw_grid % para % bo ( 1:2, 2, ip, 1 ) = get_limit (ny,rsd(2),coor(2)) 00839 pw_grid % para % bo ( 1, 3, ip, 1 ) = 1 00840 pw_grid % para % bo ( 2, 3, ip, 1 ) = nz 00841 ! distribution xYz 00842 pw_grid % para % bo ( 1:2, 1, ip, 2 ) = get_limit (nx,rsd(1),coor(1)) 00843 pw_grid % para % bo ( 1, 2, ip, 2 ) = 1 00844 pw_grid % para % bo ( 2, 2, ip, 2 ) = ny 00845 pw_grid % para % bo ( 1:2, 3, ip, 2 ) = get_limit (nz,rsd(2),coor(2)) 00846 ! distribution Xyz 00847 pw_grid % para % bo ( 1, 1, ip, 3 ) = 1 00848 pw_grid % para % bo ( 2, 1, ip, 3 ) = nx 00849 pw_grid % para % bo ( 1:2, 2, ip, 3 ) = get_limit (ny,rsd(1),coor(1)) 00850 pw_grid % para % bo ( 1:2, 3, ip, 3 ) = get_limit (nz,rsd(2),coor(2)) 00851 END DO 00852 00853 !..find the g space distribution 00854 pw_grid % ngpts_cut_local = 0 00855 00856 ALLOCATE ( pw_grid % para % nyzray ( 0: np-1 ), STAT = ierr ) 00857 CPPrecondition(ierr == 0,cp_failure_level,routineP,error,failure) 00858 00859 ALLOCATE ( pw_grid % para % yzq ( ny, nz ), STAT = ierr ) 00860 CPPrecondition(ierr == 0,cp_failure_level,routineP,error,failure) 00861 00862 IF ( pw_grid % spherical .OR. pw_grid % grid_span == HALFSPACE & 00863 .OR. .NOT. blocking ) THEN 00864 00865 pw_grid % para % ray_distribution = .TRUE. 00866 00867 pw_grid % para % yzq = -1 00868 IF ( PRESENT ( ref_grid ) ) THEN 00869 ! tag all vectors from the reference grid 00870 CALL pre_tag ( pw_grid, yz_mask, ref_grid ) 00871 END IF 00872 00873 ! Round Robin distribution 00874 ! Processors 0 .. NP-1, NP-1 .. 0 get the largest remaining batch 00875 ! of g vectors in turn 00876 00877 i1 = SIZE ( yz_mask, 1 ) 00878 i2 = SIZE ( yz_mask, 2 ) 00879 ALLOCATE ( yz_index(2,i1*i2), STAT = ierr ) 00880 CPPrecondition(ierr == 0,cp_failure_level,routineP,error,failure) 00881 CALL order_mask ( yz_mask, yz_index, error ) 00882 DO i = 1, i1*i2 00883 lo(1) = yz_index(1,i) 00884 lo(2) = yz_index(2,i) 00885 IF ( lo(1)*lo(2) == 0 ) CYCLE 00886 gmax = yz_mask ( lo(1), lo(2) ) 00887 IF ( gmax == 0 ) CYCLE 00888 yz_mask ( lo(1), lo(2) ) = 0 00889 ip = MOD ( i-1, 2*np ) 00890 IF ( ip > np - 1 ) ip = 2*np - ip - 1 00891 IF ( ip == pw_grid % para % my_pos ) THEN 00892 pw_grid % ngpts_cut_local = pw_grid % ngpts_cut_local + gmax 00893 END IF 00894 pw_grid % para % yzq ( lo(1), lo(2) ) = ip 00895 IF ( pw_grid % grid_span == HALFSPACE ) THEN 00896 m = -lo(1) - 2*lby + 2 00897 n = -lo(2) - 2*lbz + 2 00898 pw_grid % para % yzq ( m, n ) = ip 00899 yz_mask ( m, n ) = 0 00900 END IF 00901 END DO 00902 DEALLOCATE ( yz_index, STAT = ierr ) 00903 CPPrecondition(ierr == 0,cp_failure_level,routineP,error,failure) 00904 00905 ! Count the total number of rays on each processor 00906 pw_grid % para % nyzray = 0 00907 DO i = 1, nz 00908 DO j = 1, ny 00909 ip = pw_grid % para % yzq ( j, i ) 00910 IF ( ip >= 0 ) pw_grid % para % nyzray ( ip ) = & 00911 pw_grid % para % nyzray ( ip ) + 1 00912 END DO 00913 END DO 00914 00915 ! Allocate mapping array (y:z, nray, nproc) 00916 ns = MAXVAL ( pw_grid % para % nyzray ( 0: np-1 ) ) 00917 ALLOCATE ( pw_grid % para % yzp ( 2, ns, 0: np-1 ), STAT = ierr ) 00918 CPPrecondition(ierr == 0,cp_failure_level,routineP,error,failure) 00919 00920 ! Fill mapping array, recalculate nyzray for convenience 00921 pw_grid % para % nyzray = 0 00922 DO i = 1, nz 00923 DO j = 1, ny 00924 ip = pw_grid % para % yzq ( j, i ) 00925 IF ( ip >= 0 ) THEN 00926 pw_grid % para % nyzray ( ip ) = & 00927 pw_grid % para % nyzray ( ip ) + 1 00928 ns = pw_grid % para % nyzray ( ip ) 00929 pw_grid % para % yzp ( 1, ns, ip ) = j 00930 pw_grid % para % yzp ( 2, ns, ip ) = i 00931 IF ( ip == pw_grid % para % my_pos ) THEN 00932 pw_grid % para % yzq ( j, i ) = ns 00933 ELSE 00934 pw_grid % para % yzq ( j, i ) = -1 00935 END IF 00936 ELSE 00937 pw_grid % para % yzq ( j, i ) = -2 00938 END IF 00939 END DO 00940 END DO 00941 00942 pw_grid % ngpts_local = PRODUCT ( pw_grid % npts_local ) 00943 00944 ELSE 00945 ! 00946 ! block distribution of g vectors, we do not have a spherical cutoff 00947 ! 00948 00949 pw_grid % para % blocked = .TRUE. 00950 pw_grid % para % ray_distribution = .FALSE. 00951 00952 DO ip = 0, np - 1 00953 m = pw_grid % para % bo ( 2, 2, ip, 3 ) - & 00954 pw_grid % para % bo ( 1, 2, ip, 3 ) + 1 00955 n = pw_grid % para % bo ( 2, 3, ip, 3 ) - & 00956 pw_grid % para % bo ( 1, 3, ip, 3 ) + 1 00957 pw_grid % para % nyzray ( ip ) = n*m 00958 END DO 00959 00960 ipl = pw_grid % para % rs_mpo 00961 l = pw_grid % para % bo ( 2, 1, ipl, 3 ) - & 00962 pw_grid % para % bo ( 1, 1, ipl, 3 ) + 1 00963 m = pw_grid % para % bo ( 2, 2, ipl, 3 ) - & 00964 pw_grid % para % bo ( 1, 2, ipl, 3 ) + 1 00965 n = pw_grid % para % bo ( 2, 3, ipl, 3 ) - & 00966 pw_grid % para % bo ( 1, 3, ipl, 3 ) + 1 00967 pw_grid % ngpts_cut_local = l * m * n 00968 pw_grid % ngpts_local = pw_grid % ngpts_cut_local 00969 00970 pw_grid % para % yzq = 0 00971 ny = pw_grid % para % bo ( 2, 2, ipl, 3 ) - & 00972 pw_grid % para % bo ( 1, 2, ipl, 3 ) + 1 00973 DO n = pw_grid % para % bo ( 1, 3, ipl, 3 ), & 00974 pw_grid % para % bo ( 2, 3, ipl, 3 ) 00975 i = n - pw_grid % para % bo ( 1, 3, ipl, 3 ) 00976 DO m = pw_grid % para % bo ( 1, 2, ipl, 3 ), & 00977 pw_grid % para % bo ( 2, 2, ipl, 3 ) 00978 j = m - pw_grid % para % bo ( 1, 2, ipl, 3 ) + 1 00979 pw_grid % para % yzq ( m, n ) = j + i * ny 00980 END DO 00981 END DO 00982 00983 ! Allocate mapping array (y:z, nray, nproc) 00984 ns = MAXVAL ( pw_grid % para % nyzray ( 0: np-1 ) ) 00985 ALLOCATE ( pw_grid % para % yzp ( 2, ns, 0: np-1 ), STAT = ierr ) 00986 CPPrecondition(ierr == 0,cp_failure_level,routineP,error,failure) 00987 pw_grid % para % yzp = 0 00988 00989 ALLOCATE ( pemap(0:np-1), STAT = ierr ) 00990 CPPrecondition(ierr == 0,cp_failure_level,routineP,error,failure) 00991 pemap = 0 00992 pemap(pw_grid % para % my_pos) = pw_grid % para % rs_mpo 00993 CALL mp_sum(pemap,pw_grid % para % group) 00994 00995 DO ip = 0, np - 1 00996 ipp = pemap(ip) 00997 ns = 0 00998 DO n = pw_grid % para % bo ( 1, 3, ipp, 3 ), & 00999 pw_grid % para % bo ( 2, 3, ipp, 3 ) 01000 i = n - pw_grid % bounds ( 1, 3 ) + 1 01001 DO m = pw_grid % para % bo ( 1, 2, ipp, 3 ), & 01002 pw_grid % para % bo ( 2, 2, ipp, 3 ) 01003 j = m - pw_grid % bounds ( 1, 2 ) + 1 01004 ns = ns + 1 01005 pw_grid % para % yzp ( 1, ns, ip ) = j 01006 pw_grid % para % yzp ( 2, ns, ip ) = i 01007 END DO 01008 END DO 01009 CPPrecondition(ns == pw_grid%para%nyzray(ip),cp_failure_level,routineP,error,failure) 01010 END DO 01011 01012 DEALLOCATE ( pemap, STAT = ierr ) 01013 CPPrecondition(ierr == 0,cp_failure_level,routineP,error,failure) 01014 01015 END IF 01016 01017 END IF 01018 01019 ! pos_of_x(i) tells on which cpu pw%cr3d(i,:,:) is located 01020 ! should be computable in principle, without the need for communication 01021 IF (pw_grid % para % mode .EQ. PW_MODE_DISTRIBUTED) THEN 01022 ALLOCATE(pw_grid % para % pos_of_x( pw_grid % bounds(1,1): pw_grid % bounds(2,1) )) 01023 pw_grid % para % pos_of_x = 0 01024 pw_grid % para % pos_of_x(pw_grid % bounds_local(1,1) : pw_grid % bounds_local(2,1))=pw_grid % para % my_pos 01025 CALL mp_sum( pw_grid % para % pos_of_x, pw_grid % para % group ) 01026 ELSE 01027 ! this should not be needed 01028 ALLOCATE(pw_grid % para % pos_of_x( pw_grid % bounds(1,1): pw_grid % bounds(2,1) )) 01029 pw_grid % para % pos_of_x = 0 01030 ENDIF 01031 01032 CALL timestop(handle) 01033 01034 END SUBROUTINE pw_grid_distribute 01035 01036 ! ***************************************************************************** 01037 SUBROUTINE pre_tag ( pw_grid, yz_mask, ref_grid ) 01038 01039 TYPE(pw_grid_type), POINTER :: pw_grid 01040 INTEGER, DIMENSION(:, :), INTENT(INOUT) :: yz_mask 01041 TYPE(pw_grid_type), INTENT(IN) :: ref_grid 01042 01043 INTEGER :: gmax, ig, ip, lby, lbz, my, 01044 mz, ny, nz, uby, ubz, y, yp, 01045 z, zp 01046 01047 !------------------------------------------------------------------------------ 01048 01049 ny = ref_grid % npts ( 2 ) 01050 nz = ref_grid % npts ( 3 ) 01051 lby = pw_grid % bounds ( 1, 2 ) 01052 lbz = pw_grid % bounds ( 1, 3 ) 01053 uby = pw_grid % bounds ( 2, 2 ) 01054 ubz = pw_grid % bounds ( 2, 3 ) 01055 my = SIZE ( yz_mask, 1 ) 01056 mz = SIZE ( yz_mask, 2 ) 01057 01058 ! loop over all processors and all g vectors yz lines on this processor 01059 DO ip = 0, ref_grid % para % group_size - 1 01060 DO ig = 1, ref_grid % para % nyzray ( ip ) 01061 ! go from mapped coordinates to original coordinates 01062 ! 0 .. N-1 -> -n/2 .. (n+1)/2 01063 y = ref_grid % para % yzp ( 1, ig, ip ) - 1 01064 IF ( y > ny/2 ) y = y - ny 01065 z = ref_grid % para % yzp ( 2, ig, ip ) - 1 01066 IF ( z > nz/2 ) z = z - nz 01067 ! check if this is inside the realm of the new grid 01068 IF ( y < lby .OR. y > uby .OR. z < lbz .OR. z > ubz ) CYCLE 01069 ! go to shifted coordinates 01070 y = y - lby + 1 01071 z = z - lbz + 1 01072 ! this tag is outside the cutoff range of the new grid 01073 IF ( pw_grid % grid_span == HALFSPACE ) THEN 01074 yp = -y - 2*lby + 2 01075 zp = -z - 2*lbz + 2 01076 ! if the referenz grid is larger than the mirror point may be 01077 ! outside the new grid even if the original point is inside 01078 IF ( yp < 1 .OR. yp > my .OR. zp < 1 .OR. zp > mz ) CYCLE 01079 gmax = MAX ( yz_mask ( y, z ), yz_mask ( yp, zp ) ) 01080 IF ( gmax == 0 ) CYCLE 01081 yz_mask ( y, z ) = 0 01082 yz_mask ( yp, zp ) = 0 01083 pw_grid % para % yzq ( y, z ) = ip 01084 pw_grid % para % yzq ( yp, zp ) = ip 01085 ELSE 01086 gmax = yz_mask ( y, z ) 01087 IF ( gmax == 0 ) CYCLE 01088 yz_mask ( y, z ) = 0 01089 pw_grid % para % yzq ( y, z ) = ip 01090 END IF 01091 IF ( ip == pw_grid % para % my_pos ) THEN 01092 pw_grid % ngpts_cut_local = pw_grid % ngpts_cut_local + gmax 01093 END IF 01094 END DO 01095 END DO 01096 01097 END SUBROUTINE pre_tag 01098 01099 !------------------------------------------------------------------------------ 01100 01101 ! ***************************************************************************** 01102 SUBROUTINE order_mask ( yz_mask, yz_index, error ) 01103 01104 INTEGER, DIMENSION(:, :), INTENT(IN) :: yz_mask 01105 INTEGER, DIMENSION(:, :), INTENT(OUT) :: yz_index 01106 TYPE(cp_error_type), INTENT(inout) :: error 01107 01108 CHARACTER(len=*), PARAMETER :: routineN = 'order_mask', 01109 routineP = moduleN//':'//routineN 01110 01111 INTEGER :: i, i1, i2, icount, ierr, j, 01112 j1, j2 01113 INTEGER, ALLOCATABLE, DIMENSION(:) :: cindex, icol, irow, rindex 01114 LOGICAL :: failure 01115 !NB load balance 01116 INTEGER :: im, ic, jc, ii, jj 01117 01118 !------------------------------------------------------------------------------ 01119 01120 #ifdef ORIGINAL_RAY_DISTRIBUTION 01121 !NB original distribution, can have bad load balance in cp_ddapc_apply_CD which 01122 !NB uses spherical cutoff even though overall grid is full and block distributed 01123 failure = .FALSE. 01124 01125 i1 = SIZE ( yz_mask, 1 ) 01126 i2 = SIZE ( yz_mask, 2 ) 01127 ALLOCATE ( irow(i1), rindex(i1), STAT = ierr ) 01128 CPPrecondition(ierr == 0,cp_failure_level,routineP,error,failure) 01129 ALLOCATE ( icol(i2), cindex(i2), STAT = ierr ) 01130 CPPrecondition(ierr == 0,cp_failure_level,routineP,error,failure) 01131 01132 yz_index = 0 01133 DO i = 1, i1 01134 irow(i) = SUM(yz_mask(i,:)) 01135 END DO 01136 CALL sort ( irow, i1, rindex ) 01137 icount = 0 01138 DO i = i1, 1, -1 01139 j1 = rindex ( i ) 01140 icol = yz_mask(i,:) 01141 CALL sort ( icol, i2, cindex ) 01142 DO j = i2, 1, -1 01143 j2 = cindex ( j ) 01144 IF ( yz_mask(j1,j2) /= 0 ) THEN 01145 icount = icount + 1 01146 yz_index(1,icount) = j1 01147 yz_index(2,icount) = j2 01148 END IF 01149 END DO 01150 END DO 01151 01152 j=HUGE(j) 01153 DO i=1,icount 01154 j1=yz_index(1,icount) 01155 j2=yz_index(2,icount) 01156 IF(j < yz_mask(j1,j2)) STOP 01157 j = yz_mask(j1,j2) 01158 ENDDO 01159 01160 DEALLOCATE ( irow, rindex, icol, cindex, STAT = ierr ) 01161 CPPrecondition(ierr == 0,cp_failure_level,routineP,error,failure) 01162 01163 #else 01164 !NB spiral out from origin, so that even if overall grid is full and 01165 !NB block distributed, spherical cutoff still leads to good load 01166 !NB balance in cp_ddapc_apply_CD 01167 01168 i1 = SIZE ( yz_mask, 1 ) 01169 i2 = SIZE ( yz_mask, 2 ) 01170 yz_index = 0 01171 01172 icount = 1 01173 ic=i1/2 01174 jc=i2/2 01175 ii=ic 01176 jj=jc 01177 IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN 01178 IF (yz_mask(ii,jj) /= 0) THEN 01179 yz_index(1,icount) = ii 01180 yz_index(2,icount) = jj 01181 icount = icount + 1 01182 ENDIF 01183 ENDIF 01184 DO im=1, MAX(ic+1,jc+1) 01185 ii = ic-im 01186 DO jj=jc-im,jc+im 01187 IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN 01188 IF (yz_mask(ii,jj) /= 0) THEN 01189 yz_index(1,icount) = ii 01190 yz_index(2,icount) = jj 01191 icount = icount + 1 01192 ENDIF 01193 ENDIF 01194 END DO 01195 ii = ic+im 01196 DO jj=jc-im,jc+im 01197 IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN 01198 IF (yz_mask(ii,jj) /= 0) THEN 01199 yz_index(1,icount) = ii 01200 yz_index(2,icount) = jj 01201 icount = icount + 1 01202 ENDIF 01203 ENDIF 01204 END DO 01205 jj = jc-im 01206 DO ii=ic-im+1,ic+im-1 01207 IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN 01208 IF (yz_mask(ii,jj) /= 0) THEN 01209 yz_index(1,icount) = ii 01210 yz_index(2,icount) = jj 01211 icount = icount + 1 01212 ENDIF 01213 ENDIF 01214 END DO 01215 jj = jc+im 01216 DO ii=ic-im+1,ic+im-1 01217 IF (ii > 0 .AND. ii <= i1 .AND. jj > 0 .AND. jj <= i2) THEN 01218 IF (yz_mask(ii,jj) /= 0) THEN 01219 yz_index(1,icount) = ii 01220 yz_index(2,icount) = jj 01221 icount = icount + 1 01222 ENDIF 01223 ENDIF 01224 END DO 01225 END DO 01226 01227 #endif 01228 01229 END SUBROUTINE order_mask 01230 01231 ! ***************************************************************************** 01238 SUBROUTINE pw_grid_count ( h_inv, pw_grid, cutoff, yz_mask ) 01239 01240 REAL(KIND=dp), DIMENSION(3, 3) :: h_inv 01241 TYPE(pw_grid_type), POINTER :: pw_grid 01242 REAL(KIND=dp), INTENT(IN) :: cutoff 01243 INTEGER, DIMENSION(:, :), INTENT(OUT) :: yz_mask 01244 01245 CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_count', 01246 routineP = moduleN//':'//routineN 01247 01248 INTEGER :: gpt, l, m, mm, n, 01249 n_upperlimit, nlim( 2 ), nn 01250 INTEGER, DIMENSION(:, :), POINTER :: bounds 01251 REAL(KIND=dp) :: length, length_x, length_y, 01252 length_z 01253 01254 bounds => pw_grid % bounds 01255 01256 IF ( pw_grid % grid_span == HALFSPACE ) THEN 01257 n_upperlimit = 0 01258 ELSE IF ( pw_grid % grid_span == FULLSPACE ) THEN 01259 n_upperlimit = bounds ( 2, 3 ) 01260 ELSE 01261 CALL stop_program(routineN,moduleN,__LINE__,& 01262 "No type set for the grid") 01263 END IF 01264 01265 ! finds valid g-points within grid 01266 gpt = 0 01267 IF ( pw_grid % para % mode == PW_MODE_LOCAL ) THEN 01268 nlim ( 1 ) = bounds ( 1, 3 ) 01269 nlim ( 2 ) = n_upperlimit 01270 ELSE IF ( pw_grid % para % mode == PW_MODE_DISTRIBUTED ) THEN 01271 n = n_upperlimit - bounds ( 1, 3 ) + 1 01272 nlim = get_limit ( n, pw_grid % para % group_size, pw_grid % para % my_pos ) 01273 nlim = nlim + bounds ( 1, 3 ) - 1 01274 ELSE 01275 CALL stop_program(routineN,moduleN,__LINE__,& 01276 "para % mode not specified") 01277 END IF 01278 01279 yz_mask = 0 01280 DO n = nlim ( 1 ), nlim ( 2 ) 01281 nn = n - bounds ( 1, 3) + 1 01282 DO m = bounds ( 1, 2 ), bounds ( 2, 2 ) 01283 mm = m - bounds ( 1, 2) + 1 01284 DO l = bounds ( 1, 1 ), bounds ( 2, 1 ) 01285 IF ( pw_grid % grid_span == HALFSPACE .AND. n == 0 ) THEN 01286 IF ( ( m == 0 .AND. l > 0 ) .OR. ( m > 0 ) ) CYCLE 01287 END IF 01288 01289 length_x & 01290 = REAL(l,dp) * h_inv(1,1) 01291 + REAL(m,dp) * h_inv(2,1) 01292 + REAL(n,dp) * h_inv(3,1) 01293 length_y & 01294 = REAL(l,dp) * h_inv(1,2) 01295 + REAL(m,dp) * h_inv(2,2) 01296 + REAL(n,dp) * h_inv(3,2) 01297 length_z & 01298 = REAL(l,dp) * h_inv(1,3) 01299 + REAL(m,dp) * h_inv(2,3) 01300 + REAL(n,dp) * h_inv(3,3) 01301 01302 length = length_x ** 2 + length_y ** 2 + length_z ** 2 01303 length = twopi * twopi * length 01304 01305 IF ( 0.5_dp * length <= cutoff ) THEN 01306 gpt = gpt + 1 01307 yz_mask ( mm, nn ) = yz_mask ( mm, nn ) + 1 01308 END IF 01309 01310 END DO 01311 END DO 01312 END DO 01313 01314 ! number of g-vectors for grid 01315 IF ( pw_grid % para % mode == PW_MODE_DISTRIBUTED ) THEN 01316 CALL mp_sum ( gpt, pw_grid % para % group ) 01317 CALL mp_sum ( yz_mask, pw_grid % para % group ) 01318 ENDIF 01319 pw_grid % ngpts_cut = gpt 01320 01321 END SUBROUTINE pw_grid_count 01322 01323 ! ***************************************************************************** 01330 SUBROUTINE pw_grid_assign ( h_inv, pw_grid, cutoff, error ) 01331 01332 REAL(KIND=dp), DIMENSION(3, 3) :: h_inv 01333 TYPE(pw_grid_type), POINTER :: pw_grid 01334 REAL(KIND=dp), INTENT(IN) :: cutoff 01335 TYPE(cp_error_type), INTENT(inout) :: error 01336 01337 CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_assign', 01338 routineP = moduleN//':'//routineN 01339 01340 INTEGER :: gpt, handle, i, ip, l, lby, 01341 lbz, ll, m, mm, n, 01342 n_upperlimit, nn 01343 INTEGER, DIMENSION(2, 3) :: bol, bounds 01344 LOGICAL :: failure 01345 REAL(KIND=dp) :: length, length_x, length_y, 01346 length_z 01347 01348 CALL timeset(routineN,handle) 01349 01350 failure = .FALSE. 01351 01352 bounds = pw_grid % bounds 01353 lby = pw_grid % bounds ( 1, 2 ) 01354 lbz = pw_grid % bounds ( 1, 3 ) 01355 01356 IF ( pw_grid % grid_span == HALFSPACE ) THEN 01357 n_upperlimit = 0 01358 ELSE IF ( pw_grid % grid_span == FULLSPACE ) THEN 01359 n_upperlimit = bounds ( 2, 3 ) 01360 ELSE 01361 CALL stop_program(routineN,moduleN,__LINE__,& 01362 "No type set for the grid") 01363 END IF 01364 01365 ! finds valid g-points within grid 01366 IF ( pw_grid % para % mode == PW_MODE_LOCAL ) THEN 01367 gpt = 0 01368 DO n = bounds ( 1, 3 ), n_upperlimit 01369 DO m = bounds ( 1, 2 ), bounds ( 2, 2 ) 01370 DO l = bounds ( 1, 1 ), bounds ( 2, 1 ) 01371 IF ( pw_grid % grid_span == HALFSPACE .AND. n == 0 ) THEN 01372 IF ( ( m == 0 .AND. l > 0 ) .OR. ( m > 0 ) ) CYCLE 01373 END IF 01374 01375 length_x & 01376 = REAL(l,dp) * h_inv(1,1) 01377 + REAL(m,dp) * h_inv(2,1) 01378 + REAL(n,dp) * h_inv(3,1) 01379 length_y & 01380 = REAL(l,dp) * h_inv(1,2) 01381 + REAL(m,dp) * h_inv(2,2) 01382 + REAL(n,dp) * h_inv(3,2) 01383 length_z & 01384 = REAL(l,dp) * h_inv(1,3) 01385 + REAL(m,dp) * h_inv(2,3) 01386 + REAL(n,dp) * h_inv(3,3) 01387 01388 length = length_x ** 2 + length_y ** 2 + length_z ** 2 01389 length = twopi * twopi * length 01390 01391 IF ( 0.5_dp * length <= cutoff ) THEN 01392 gpt = gpt + 1 01393 pw_grid % g ( 1, gpt ) = twopi * length_x 01394 pw_grid % g ( 2, gpt ) = twopi * length_y 01395 pw_grid % g ( 3, gpt ) = twopi * length_z 01396 pw_grid % gsq ( gpt ) = length 01397 pw_grid % g_hat ( 1, gpt ) = l 01398 pw_grid % g_hat ( 2, gpt ) = m 01399 pw_grid % g_hat ( 3, gpt ) = n 01400 END IF 01401 01402 END DO 01403 END DO 01404 END DO 01405 01406 ELSE 01407 01408 IF ( pw_grid % para % ray_distribution ) THEN 01409 01410 gpt = 0 01411 ip = pw_grid % para % my_pos 01412 DO i = 1, pw_grid % para % nyzray ( ip ) 01413 n = pw_grid % para % yzp ( 2, i, ip ) + lbz - 1 01414 m = pw_grid % para % yzp ( 1, i, ip ) + lby - 1 01415 IF ( n > n_upperlimit ) CYCLE 01416 DO l = bounds ( 1, 1 ), bounds ( 2, 1 ) 01417 IF ( pw_grid % grid_span == HALFSPACE .AND. n == 0 ) THEN 01418 IF ( ( m == 0 .AND. l > 0 ) .OR. ( m > 0 ) ) CYCLE 01419 END IF 01420 01421 length_x & 01422 = REAL(l,dp) * h_inv(1,1) 01423 + REAL(m,dp) * h_inv(2,1) 01424 + REAL(n,dp) * h_inv(3,1) 01425 length_y & 01426 = REAL(l,dp) * h_inv(1,2) 01427 + REAL(m,dp) * h_inv(2,2) 01428 + REAL(n,dp) * h_inv(3,2) 01429 length_z & 01430 = REAL(l,dp) * h_inv(1,3) 01431 + REAL(m,dp) * h_inv(2,3) 01432 + REAL(n,dp) * h_inv(3,3) 01433 01434 length = length_x ** 2 + length_y ** 2 + length_z ** 2 01435 length = twopi * twopi * length 01436 01437 IF ( 0.5_dp * length <= cutoff ) THEN 01438 gpt = gpt + 1 01439 pw_grid % g ( 1, gpt ) = twopi * length_x 01440 pw_grid % g ( 2, gpt ) = twopi * length_y 01441 pw_grid % g ( 3, gpt ) = twopi * length_z 01442 pw_grid % gsq ( gpt ) = length 01443 pw_grid % g_hat ( 1, gpt ) = l 01444 pw_grid % g_hat ( 2, gpt ) = m 01445 pw_grid % g_hat ( 3, gpt ) = n 01446 END IF 01447 01448 END DO 01449 END DO 01450 01451 ELSE 01452 01453 bol = pw_grid % para % bo ( :, :, pw_grid % para % rs_mpo, 3 ) 01454 gpt = 0 01455 DO n = bounds ( 1, 3 ), bounds ( 2, 3 ) 01456 IF ( n < 0 ) THEN 01457 nn = n + pw_grid % npts ( 3 ) + 1 01458 ELSE 01459 nn = n + 1 01460 END IF 01461 IF ( nn < bol ( 1, 3 ) .OR. nn > bol ( 2, 3 ) ) CYCLE 01462 DO m = bounds ( 1, 2 ), bounds ( 2, 2 ) 01463 IF ( m < 0 ) THEN 01464 mm = m + pw_grid % npts ( 2 ) + 1 01465 ELSE 01466 mm = m + 1 01467 END IF 01468 IF ( mm < bol ( 1, 2 ) .OR. mm > bol ( 2, 2 ) ) CYCLE 01469 DO l = bounds ( 1, 1 ), bounds ( 2, 1 ) 01470 IF ( l < 0 ) THEN 01471 ll = l + pw_grid % npts ( 1 ) + 1 01472 ELSE 01473 ll = l + 1 01474 END IF 01475 IF ( ll < bol ( 1, 1 ) .OR. ll > bol ( 2, 1 ) ) CYCLE 01476 01477 length_x & 01478 = REAL(l,dp) * h_inv(1,1) 01479 + REAL(m,dp) * h_inv(2,1) 01480 + REAL(n,dp) * h_inv(3,1) 01481 length_y & 01482 = REAL(l,dp) * h_inv(1,2) 01483 + REAL(m,dp) * h_inv(2,2) 01484 + REAL(n,dp) * h_inv(3,2) 01485 length_z & 01486 = REAL(l,dp) * h_inv(1,3) 01487 + REAL(m,dp) * h_inv(2,3) 01488 + REAL(n,dp) * h_inv(3,3) 01489 01490 length = length_x ** 2 + length_y ** 2 + length_z ** 2 01491 length = twopi * twopi * length 01492 01493 gpt = gpt + 1 01494 pw_grid % g ( 1, gpt ) = twopi * length_x 01495 pw_grid % g ( 2, gpt ) = twopi * length_y 01496 pw_grid % g ( 3, gpt ) = twopi * length_z 01497 pw_grid % gsq ( gpt ) = length 01498 pw_grid % g_hat ( 1, gpt ) = l 01499 pw_grid % g_hat ( 2, gpt ) = m 01500 pw_grid % g_hat ( 3, gpt ) = n 01501 01502 END DO 01503 END DO 01504 END DO 01505 01506 END IF 01507 01508 END IF 01509 01510 ! Check the number of g-vectors for grid 01511 CPPrecondition(pw_grid % ngpts_cut_local == gpt,cp_failure_level,routineP,error,failure) 01512 IF ( pw_grid % para % mode == PW_MODE_DISTRIBUTED ) THEN 01513 CALL mp_sum ( gpt, pw_grid % para % group ) 01514 CPPrecondition(pw_grid % ngpts_cut == gpt,cp_failure_level,routineP,error,failure) 01515 ENDIF 01516 01517 pw_grid % have_g0 = .FALSE. 01518 pw_grid % first_gne0 = 1 01519 DO gpt = 1, pw_grid % ngpts_cut_local 01520 IF ( ALL ( pw_grid % g_hat ( :, gpt ) == 0 ) ) THEN 01521 pw_grid % have_g0 = .TRUE. 01522 pw_grid % first_gne0 = 2 01523 EXIT 01524 END IF 01525 END DO 01526 01527 CALL pw_grid_set_maps ( pw_grid % grid_span, pw_grid % g_hat, & 01528 pw_grid % mapl, pw_grid % mapm, pw_grid % mapn, pw_grid % npts ) 01529 01530 CALL timestop(handle) 01531 01532 END SUBROUTINE pw_grid_assign 01533 01534 ! ***************************************************************************** 01543 SUBROUTINE pw_grid_set_maps ( grid_span, g_hat, mapl, mapm, mapn, npts ) 01544 01545 INTEGER, INTENT(IN) :: grid_span 01546 INTEGER, DIMENSION(:, :), INTENT(IN) :: g_hat 01547 TYPE(map_pn), INTENT(INOUT) :: mapl, mapm, mapn 01548 INTEGER, DIMENSION(:), INTENT(IN) :: npts 01549 01550 INTEGER :: gpt, l, m, n, ng 01551 01552 ng = SIZE ( g_hat ,2 ) 01553 01554 DO gpt = 1, ng 01555 l = g_hat ( 1, gpt ) 01556 m = g_hat ( 2, gpt ) 01557 n = g_hat ( 3, gpt ) 01558 IF ( l < 0 ) THEN 01559 mapl % pos ( l ) = l + npts ( 1 ) 01560 ELSE 01561 mapl % pos ( l ) = l 01562 END IF 01563 IF ( m < 0 ) THEN 01564 mapm % pos ( m ) = m + npts ( 2 ) 01565 ELSE 01566 mapm % pos ( m ) = m 01567 END IF 01568 IF ( n < 0 ) THEN 01569 mapn % pos ( n ) = n + npts ( 3 ) 01570 ELSE 01571 mapn % pos ( n ) = n 01572 END IF 01573 01574 ! Generating the maps to the full 3-d space 01575 01576 IF ( grid_span == HALFSPACE ) THEN 01577 01578 IF ( l <= 0 ) THEN 01579 mapl % neg ( l ) = - l 01580 ELSE 01581 mapl % neg ( l ) = npts ( 1 ) - l 01582 END IF 01583 IF ( m <= 0 ) THEN 01584 mapm % neg ( m ) = - m 01585 ELSE 01586 mapm % neg ( m ) = npts ( 2 ) - m 01587 END IF 01588 IF ( n <= 0 ) THEN 01589 mapn % neg ( n ) = - n 01590 ELSE 01591 mapn % neg ( n ) = npts ( 3 ) - n 01592 END IF 01593 01594 END IF 01595 01596 END DO 01597 01598 END SUBROUTINE pw_grid_set_maps 01599 01600 ! ***************************************************************************** 01609 SUBROUTINE pw_grid_allocate ( pw_grid, ng, bounds, error ) 01610 01611 ! Argument 01612 TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid 01613 INTEGER, INTENT(IN) :: ng 01614 INTEGER, DIMENSION(:, :), INTENT(IN) :: bounds 01615 TYPE(cp_error_type), INTENT(inout) :: error 01616 01617 CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_allocate', 01618 routineP = moduleN//':'//routineN 01619 01620 INTEGER :: allocstat, handle 01621 LOGICAL :: failure 01622 01623 CALL timeset(routineN,handle) 01624 01625 failure = .FALSE. 01626 01627 ALLOCATE ( pw_grid % g ( 3, ng ), STAT = allocstat ) 01628 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 01629 ALLOCATE ( pw_grid % gsq ( ng ), STAT = allocstat ) 01630 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 01631 ALLOCATE ( pw_grid % g_hat ( 3, ng ), STAT = allocstat ) 01632 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 01633 01634 IF (pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN 01635 ALLOCATE ( pw_grid % grays ( pw_grid%npts(1), & 01636 pw_grid%para%nyzray(pw_grid%para%my_pos) ), & 01637 STAT = allocstat ) 01638 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 01639 END IF 01640 01641 ALLOCATE ( pw_grid % mapl % pos ( bounds ( 1, 1 ):bounds ( 2, 1 ) ), & 01642 STAT = allocstat ) 01643 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 01644 ALLOCATE ( pw_grid % mapl % neg ( bounds ( 1, 1 ):bounds ( 2, 1 ) ), & 01645 STAT = allocstat ) 01646 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 01647 01648 ALLOCATE ( pw_grid % mapm % pos ( bounds ( 1, 2 ):bounds ( 2, 2 ) ), & 01649 STAT = allocstat ) 01650 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 01651 ALLOCATE ( pw_grid % mapm % neg ( bounds ( 1, 2 ):bounds ( 2, 2 ) ), & 01652 STAT = allocstat ) 01653 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 01654 01655 ALLOCATE ( pw_grid % mapn % pos ( bounds ( 1, 3 ):bounds ( 2, 3 ) ), & 01656 STAT = allocstat ) 01657 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 01658 ALLOCATE ( pw_grid % mapn % neg ( bounds ( 1, 3 ):bounds ( 2, 3 ) ), & 01659 STAT = allocstat ) 01660 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 01661 01662 CALL timestop(handle) 01663 01664 END SUBROUTINE pw_grid_allocate 01665 01666 ! ***************************************************************************** 01683 SUBROUTINE pw_grid_sort ( pw_grid, ref_grid, error ) 01684 01685 ! Argument 01686 TYPE(pw_grid_type), INTENT(INOUT) :: pw_grid 01687 TYPE(pw_grid_type), INTENT(IN), OPTIONAL :: ref_grid 01688 TYPE(cp_error_type), INTENT(inout) :: error 01689 01690 CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_sort', 01691 routineP = moduleN//':'//routineN 01692 01693 INTEGER :: allocstat, handle, i, ig, ih, 01694 ip, is, it, ng, ngr 01695 INTEGER, ALLOCATABLE, DIMENSION(:) :: idx 01696 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: int_tmp 01697 LOGICAL :: failure, g_found 01698 REAL(KIND=dp) :: gig, gigr 01699 REAL(KIND=dp), ALLOCATABLE, 01700 DIMENSION(:, :) :: real_tmp 01701 01702 CALL timeset(routineN,handle) 01703 01704 failure = .FALSE. 01705 01706 ng = SIZE ( pw_grid % gsq ) 01707 ALLOCATE ( idx ( ng ), STAT = allocstat ) 01708 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 01709 01710 ! grids are (locally) ordered by length of G-vectors 01711 CALL sort ( pw_grid % gsq, ng, idx ) 01712 ! within shells order wrt x,y,z 01713 CALL sort_shells ( pw_grid % gsq, pw_grid % g_hat, idx, error ) 01714 01715 ALLOCATE ( real_tmp ( 3, ng ), STAT = allocstat ) 01716 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 01717 DO i=1,ng 01718 real_tmp(1,i)= pw_grid % g (1,idx (i)) 01719 real_tmp(2,i)= pw_grid % g (2,idx (i)) 01720 real_tmp(3,i)= pw_grid % g (3,idx (i)) 01721 ENDDO 01722 DO i=1,ng 01723 pw_grid % g (1,i)=real_tmp(1,i) 01724 pw_grid % g (2,i)=real_tmp(2,i) 01725 pw_grid % g (3,i)=real_tmp(3,i) 01726 ENDDO 01727 DEALLOCATE ( real_tmp, STAT = allocstat ) 01728 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 01729 01730 ALLOCATE ( int_tmp ( 3, ng ), STAT = allocstat ) 01731 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 01732 DO i=1,ng 01733 int_tmp(1,i)= pw_grid % g_hat (1,idx (i)) 01734 int_tmp(2,i)= pw_grid % g_hat (2,idx (i)) 01735 int_tmp(3,i)= pw_grid % g_hat (3,idx (i)) 01736 ENDDO 01737 DO i=1,ng 01738 pw_grid % g_hat (1,i)=int_tmp(1,i) 01739 pw_grid % g_hat (2,i)=int_tmp(2,i) 01740 pw_grid % g_hat (3,i)=int_tmp(3,i) 01741 ENDDO 01742 DEALLOCATE ( int_tmp, STAT = allocstat ) 01743 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 01744 01745 DEALLOCATE ( idx, STAT = allocstat ) 01746 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 01747 01748 ! check if ordering is compatible to reference grid 01749 IF ( PRESENT ( ref_grid ) ) THEN 01750 ngr = SIZE ( ref_grid % gsq ) 01751 ngr = MIN ( ng, ngr ) 01752 IF ( pw_grid % spherical ) THEN 01753 IF ( .NOT. ALL ( pw_grid % g_hat ( 1:3, 1:ngr ) & 01754 == ref_grid % g_hat ( 1:3, 1:ngr ) ) ) THEN 01755 CALL stop_program(routineN,moduleN,__LINE__,& 01756 "G space sorting not compatible") 01757 END IF 01758 ELSE 01759 ALLOCATE ( pw_grid%gidx(1:ngr), STAT = allocstat ) 01760 CPPrecondition(allocstat == 0,cp_failure_level,routineP,error,failure) 01761 pw_grid%gidx = 0 01762 ! first try as many trivial associations as possible 01763 it = 0 01764 DO ig = 1, ngr 01765 IF ( .NOT. ALL ( pw_grid % g_hat ( 1:3, ig ) & 01766 == ref_grid % g_hat ( 1:3, ig ) ) ) EXIT 01767 pw_grid%gidx(ig) = ig 01768 it = ig 01769 END DO 01770 ! now for the ones that could not be done 01771 IF ( ng == ngr ) THEN 01772 ! for the case pw_grid <= ref_grid 01773 is = it 01774 DO ig = it + 1, ngr 01775 gig = pw_grid % gsq(ig) 01776 gigr = MAX(1._dp,gig) 01777 g_found=.FALSE. 01778 DO ih = is + 1, SIZE ( ref_grid % gsq ) 01779 IF ( ABS ( gig - ref_grid % gsq(ih) )/gigr > 1.e-12_dp ) CYCLE 01780 g_found=.TRUE. 01781 EXIT 01782 END DO 01783 IF ( .NOT. g_found ) THEN 01784 WRITE(*,"(A,I10,F20.10)") "G-vector", ig, pw_grid % gsq(ig) 01785 CALL stop_program(routineN,moduleN,__LINE__,& 01786 "G vector not found") 01787 END IF 01788 ip = ih - 1 01789 DO ih = ip + 1, SIZE ( ref_grid % gsq ) 01790 IF ( ABS ( gig - ref_grid % gsq(ih) )/gigr > 1.e-12_dp ) CYCLE 01791 IF ( pw_grid % g_hat(1,ig) /= ref_grid % g_hat(1,ih) ) CYCLE 01792 IF ( pw_grid % g_hat(2,ig) /= ref_grid % g_hat(2,ih) ) CYCLE 01793 IF ( pw_grid % g_hat(3,ig) /= ref_grid % g_hat(3,ih) ) CYCLE 01794 pw_grid%gidx(ig) = ih 01795 EXIT 01796 END DO 01797 IF ( pw_grid%gidx(ig) == 0 ) THEN 01798 WRITE ( *, "(A,2I10)" ) " G-Shell ",is+1,ip+1 01799 WRITE ( *, "(A,I10,3I6,F20.10)" ) & 01800 " G-vector", ig, pw_grid % g_hat(1:3,ig), pw_grid % gsq(ig) 01801 DO ih = 1, SIZE ( ref_grid % gsq ) 01802 IF ( pw_grid % g_hat(1,ig) /= ref_grid % g_hat(1,ih) ) CYCLE 01803 IF ( pw_grid % g_hat(2,ig) /= ref_grid % g_hat(2,ih) ) CYCLE 01804 IF ( pw_grid % g_hat(3,ig) /= ref_grid % g_hat(3,ih) ) CYCLE 01805 WRITE ( *, "(A,I10,3I6,F20.10)" ) & 01806 " G-vector", ih, ref_grid % g_hat(1:3,ih),ref_grid % gsq(ih) 01807 END DO 01808 CALL stop_program(routineN,moduleN,__LINE__,"G vector not found") 01809 END IF 01810 is = ip 01811 END DO 01812 ELSE 01813 ! for the case pw_grid > ref_grid 01814 is = it 01815 DO ig = it + 1, ngr 01816 gig = ref_grid % gsq(ig) 01817 gigr = MAX(1._dp,gig) 01818 g_found=.FALSE. 01819 DO ih = is + 1, ng 01820 IF ( ABS ( pw_grid % gsq(ih) - gig )/gigr > 1.e-12_dp ) CYCLE 01821 g_found=.TRUE. 01822 EXIT 01823 END DO 01824 IF ( .NOT. g_found ) THEN 01825 WRITE(*,"(A,I10,F20.10)") "G-vector", ig, ref_grid % gsq(ig) 01826 CALL stop_program(routineN,moduleN,__LINE__,"G vector not found") 01827 END IF 01828 ip = ih - 1 01829 DO ih = ip + 1, ng 01830 IF ( ABS ( pw_grid % gsq(ih) - gig )/gigr > 1.e-12_dp ) CYCLE 01831 IF ( pw_grid % g_hat(1,ih) /= ref_grid % g_hat(1,ig) ) CYCLE 01832 IF ( pw_grid % g_hat(2,ih) /= ref_grid % g_hat(2,ig) ) CYCLE 01833 IF ( pw_grid % g_hat(3,ih) /= ref_grid % g_hat(3,ig) ) CYCLE 01834 pw_grid%gidx(ig) = ih 01835 EXIT 01836 END DO 01837 IF ( pw_grid%gidx(ig) == 0 ) THEN 01838 WRITE ( *, "(A,2I10)" ) " G-Shell ",is+1,ip+1 01839 WRITE ( *, "(A,I10,3I6,F20.10)" ) & 01840 " G-vector", ig, ref_grid % g_hat(1:3,ig), ref_grid % gsq(ig) 01841 DO ih = 1, ng 01842 IF ( pw_grid % g_hat(1,ih) /= ref_grid % g_hat(1,ig) ) CYCLE 01843 IF ( pw_grid % g_hat(2,ih) /= ref_grid % g_hat(2,ig) ) CYCLE 01844 IF ( pw_grid % g_hat(3,ih) /= ref_grid % g_hat(3,ig) ) CYCLE 01845 WRITE ( *, "(A,I10,3I6,F20.10)" ) & 01846 " G-vector", ih, pw_grid % g_hat(1:3,ih),pw_grid % gsq(ih) 01847 END DO 01848 CALL stop_program(routineN,moduleN,__LINE__,"G vector not found") 01849 END IF 01850 is = ip 01851 END DO 01852 END IF 01853 ! test if all g-vectors are associated 01854 IF ( ANY ( pw_grid%gidx == 0 ) ) THEN 01855 CALL stop_program(routineN,moduleN,__LINE__,"G space sorting not compatible") 01856 END IF 01857 END IF 01858 END IF 01859 01860 !check if G=0 is at first position 01861 IF ( pw_grid % have_g0 ) THEN 01862 IF ( pw_grid % gsq ( 1 ) /= 0.0_dp ) THEN 01863 CALL stop_program(routineN,moduleN,__LINE__,"G = 0 not in first position") 01864 END IF 01865 END IF 01866 01867 CALL timestop(handle) 01868 01869 END SUBROUTINE pw_grid_sort 01870 01871 ! ***************************************************************************** 01872 SUBROUTINE sort_shells ( gsq, g_hat, idx, error ) 01873 01874 ! Argument 01875 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: gsq 01876 INTEGER, DIMENSION(:, :), INTENT(IN) :: g_hat 01877 INTEGER, DIMENSION(:), INTENT(INOUT) :: idx 01878 TYPE(cp_error_type), INTENT(inout) :: error 01879 01880 CHARACTER(len=*), PARAMETER :: routineN = 'sort_shells', 01881 routineP = moduleN//':'//routineN 01882 REAL(KIND=dp), PARAMETER :: small = 5.e-16_dp 01883 01884 INTEGER :: handle, ig, ng, s1, s2 01885 LOGICAL :: failure 01886 REAL(KIND=dp) :: s_begin 01887 01888 CALL timeset(routineN,handle) 01889 01890 ! Juergs temporary hack to get the grids sorted for large (4000Ry) cutoffs. 01891 ! might need to call lapack for machine precision. 01892 01893 failure = .FALSE. 01894 01895 ng = SIZE ( gsq ) 01896 s_begin = -1.0_dp 01897 s1 = 0 01898 s2 = 0 01899 ig = 1 01900 DO ig = 1, ng 01901 IF ( ABS ( gsq ( ig ) - s_begin ) < small ) THEN 01902 s2 = ig 01903 ELSE 01904 CALL redist ( g_hat, idx, s1, s2, error) 01905 s_begin = gsq ( ig ) 01906 s1 = ig 01907 s2 = ig 01908 END IF 01909 END DO 01910 CALL redist ( g_hat, idx, s1, s2, error ) 01911 01912 CALL timestop(handle) 01913 01914 END SUBROUTINE sort_shells 01915 01916 ! ***************************************************************************** 01917 SUBROUTINE redist ( g_hat, idx, s1, s2, error ) 01918 01919 ! Argument 01920 INTEGER, DIMENSION(:, :), INTENT(IN) :: g_hat 01921 INTEGER, DIMENSION(:), INTENT(INOUT) :: idx 01922 INTEGER, INTENT(IN) :: s1, s2 01923 TYPE(cp_error_type), INTENT(inout) :: error 01924 01925 CHARACTER(len=*), PARAMETER :: routineN = 'redist', 01926 routineP = moduleN//':'//routineN 01927 01928 INTEGER :: i, ii, info, n1, n2, n3, ns 01929 INTEGER, ALLOCATABLE, DIMENSION(:) :: indl 01930 LOGICAL :: failure 01931 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: slen 01932 01933 failure = .FALSE. 01934 01935 IF ( s2 <= s1 ) RETURN 01936 ns = s2 - s1 + 1 01937 ALLOCATE ( indl ( ns ), STAT = info ) 01938 CPPrecondition(info == 0,cp_failure_level,routineP,error,failure) 01939 ALLOCATE ( slen ( ns ), STAT = info ) 01940 CPPrecondition(info == 0,cp_failure_level,routineP,error,failure) 01941 01942 DO i = s1, s2 01943 ii = idx ( i ) 01944 n1 = g_hat(1,ii) 01945 n2 = g_hat(2,ii) 01946 n3 = g_hat(3,ii) 01947 slen ( i - s1 + 1 ) = 1000.0_dp * REAL(n1,dp) + 01948 REAL(n2,dp) + 0.001_dp * REAL(n3,dp) 01949 END DO 01950 CALL sort ( slen, ns, indl ) 01951 DO i = 1, ns 01952 ii = indl ( i ) + s1 - 1 01953 indl ( i ) = idx ( ii ) 01954 END DO 01955 idx ( s1:s2 ) = indl ( 1:ns ) 01956 01957 DEALLOCATE ( indl, STAT = info ) 01958 CPPrecondition(info == 0,cp_failure_level,routineP,error,failure) 01959 DEALLOCATE ( slen, STAT = info ) 01960 CPPrecondition(info == 0,cp_failure_level,routineP,error,failure) 01961 01962 END SUBROUTINE redist 01963 01964 ! ***************************************************************************** 01970 SUBROUTINE pw_grid_remap ( pw_grid, yz, error ) 01971 01972 ! Argument 01973 TYPE(pw_grid_type), POINTER :: pw_grid 01974 INTEGER, DIMENSION(:, :), INTENT(OUT) :: yz 01975 TYPE(cp_error_type), INTENT(inout) :: error 01976 01977 CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_remap', 01978 routineP = moduleN//':'//routineN 01979 01980 INTEGER :: gpt, handle, i, ip, is, j, m, 01981 n, ny, nz 01982 INTEGER, DIMENSION(:), POINTER :: mapm, mapn 01983 LOGICAL :: failure 01984 01985 failure = .FALSE. 01986 01987 IF ( pw_grid % para % mode == PW_MODE_LOCAL ) RETURN 01988 01989 CALL timeset(routineN,handle) 01990 01991 ny = pw_grid % npts ( 2 ) 01992 nz = pw_grid % npts ( 3 ) 01993 01994 yz = 0 01995 pw_grid % para % yzp = 0 01996 pw_grid % para % yzq = 0 01997 01998 mapm => pw_grid % mapm % pos 01999 mapn => pw_grid % mapn % pos 02000 02001 DO gpt = 1, SIZE ( pw_grid % gsq ) 02002 m = mapm ( pw_grid % g_hat ( 2, gpt ) ) + 1 02003 n = mapn ( pw_grid % g_hat ( 3, gpt ) ) + 1 02004 yz ( m, n ) = yz ( m, n ) + 1 02005 END DO 02006 IF ( pw_grid % grid_span == HALFSPACE ) THEN 02007 mapm => pw_grid % mapm % neg 02008 mapn => pw_grid % mapn % neg 02009 DO gpt = 1, SIZE ( pw_grid % gsq ) 02010 m = mapm ( pw_grid % g_hat ( 2, gpt ) ) + 1 02011 n = mapn ( pw_grid % g_hat ( 3, gpt ) ) + 1 02012 yz ( m, n ) = yz ( m, n ) + 1 02013 END DO 02014 END IF 02015 02016 ip = pw_grid % para % my_pos 02017 is = 0 02018 DO i = 1, nz 02019 DO j = 1, ny 02020 IF ( yz ( j, i ) > 0 ) THEN 02021 is = is + 1 02022 pw_grid % para % yzp ( 1, is, ip ) = j 02023 pw_grid % para % yzp ( 2, is, ip ) = i 02024 pw_grid % para % yzq ( j, i ) = is 02025 END IF 02026 END DO 02027 END DO 02028 02029 CPPrecondition(is == pw_grid % para % nyzray ( ip ),cp_failure_level,routineP,error,failure) 02030 CALL mp_sum ( pw_grid % para % yzp, pw_grid % para % group ) 02031 02032 CALL timestop(handle) 02033 02034 END SUBROUTINE pw_grid_remap 02035 02036 ! ***************************************************************************** 02047 SUBROUTINE pw_grid_change ( cell, pw_grid ) 02048 ! Argument 02049 TYPE(cell_type), POINTER :: cell 02050 TYPE(pw_grid_type), POINTER :: pw_grid 02051 02052 INTEGER :: gpt 02053 REAL(KIND=dp) :: l, m, n 02054 REAL(KIND=dp), DIMENSION(:, :), POINTER :: g, h_inv 02055 02056 h_inv => cell % h_inv 02057 g => pw_grid % g 02058 02059 pw_grid % vol = ABS ( cell % deth ) 02060 pw_grid % dvol = pw_grid % vol / REAL ( pw_grid % ngpts,KIND=dp) 02061 pw_grid % dr ( 1 ) = SQRT ( SUM ( cell % hmat ( :, 1 ) ** 2 ) ) & 02062 / REAL ( pw_grid % npts ( 1 ),KIND=dp) 02063 pw_grid % dr ( 2 ) = SQRT ( SUM ( cell % hmat ( :, 2 ) ** 2 ) ) & 02064 / REAL ( pw_grid % npts ( 2 ),KIND=dp) 02065 pw_grid % dr ( 3 ) = SQRT ( SUM ( cell % hmat ( :, 3 ) ** 2 ) ) & 02066 / REAL ( pw_grid % npts ( 3 ),KIND=dp) 02067 pw_grid % dh ( :,1 ) = cell % hmat ( :, 1 ) / REAL ( pw_grid % npts ( 1 ),KIND=dp) 02068 pw_grid % dh ( :,2 ) = cell % hmat ( :, 2 ) / REAL ( pw_grid % npts ( 2 ),KIND=dp) 02069 pw_grid % dh ( :,3 ) = cell % hmat ( :, 3 ) / REAL ( pw_grid % npts ( 3 ),KIND=dp) 02070 pw_grid % dh_inv ( 1,: ) = cell % h_inv ( 1,: ) * REAL ( pw_grid % npts ( 1 ),KIND=dp) 02071 pw_grid % dh_inv ( 2,: ) = cell % h_inv ( 2,: ) * REAL ( pw_grid % npts ( 2 ),KIND=dp) 02072 pw_grid % dh_inv ( 3,: ) = cell % h_inv ( 3,: ) * REAL ( pw_grid % npts ( 3 ),KIND=dp) 02073 pw_grid % orthorhombic = cell % orthorhombic 02074 02075 DO gpt = 1, SIZE ( g ,2 ) 02076 02077 l = twopi * REAL ( pw_grid % g_hat ( 1, gpt ),KIND=dp) 02078 m = twopi * REAL ( pw_grid % g_hat ( 2, gpt ),KIND=dp) 02079 n = twopi * REAL ( pw_grid % g_hat ( 3, gpt ),KIND=dp) 02080 02081 g ( 1, gpt ) = l * h_inv(1,1) + m * h_inv(2,1) + n * h_inv(3,1) 02082 g ( 2, gpt ) = l * h_inv(1,2) + m * h_inv(2,2) + n * h_inv(3,2) 02083 g ( 3, gpt ) = l * h_inv(1,3) + m * h_inv(2,3) + n * h_inv(3,3) 02084 02085 pw_grid % gsq ( gpt ) = g ( 1, gpt ) * g ( 1, gpt ) & 02086 + g ( 2, gpt ) * g ( 2, gpt ) & 02087 + g ( 3, gpt ) * g ( 3, gpt ) 02088 02089 END DO 02090 02091 END SUBROUTINE pw_grid_change 02092 02093 ! ***************************************************************************** 02104 SUBROUTINE pw_grid_retain(pw_grid, error) 02105 TYPE(pw_grid_type), POINTER :: pw_grid 02106 TYPE(cp_error_type), INTENT(inout) :: error 02107 02108 CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_retain', 02109 routineP = moduleN//':'//routineN 02110 02111 LOGICAL :: failure 02112 02113 failure=.FALSE. 02114 02115 CPPrecondition(ASSOCIATED(pw_grid),cp_failure_level,routineP,error,failure) 02116 IF (.NOT. failure) THEN 02117 CPPreconditionNoFail(pw_grid%ref_count>0,cp_failure_level,routineP,error) 02118 pw_grid%ref_count=pw_grid%ref_count+1 02119 END IF 02120 END SUBROUTINE pw_grid_retain 02121 02122 ! ***************************************************************************** 02133 SUBROUTINE pw_grid_release(pw_grid, error) 02134 TYPE(pw_grid_type), POINTER :: pw_grid 02135 TYPE(cp_error_type), INTENT(inout) :: error 02136 02137 CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_release', 02138 routineP = moduleN//':'//routineN 02139 02140 INTEGER :: stat 02141 LOGICAL :: failure 02142 02143 failure=.FALSE. 02144 02145 IF (ASSOCIATED(pw_grid)) THEN 02146 CPPreconditionNoFail(pw_grid%ref_count>0,cp_failure_level,routineP,error) 02147 pw_grid%ref_count=pw_grid%ref_count-1 02148 IF (pw_grid%ref_count==0) THEN 02149 IF ( ASSOCIATED ( pw_grid % gidx ) ) THEN 02150 DEALLOCATE ( pw_grid % gidx, STAT = stat ) 02151 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02152 END IF 02153 IF ( ASSOCIATED ( pw_grid % g ) ) THEN 02154 DEALLOCATE ( pw_grid % g, STAT = stat ) 02155 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02156 END IF 02157 IF ( ASSOCIATED ( pw_grid % gsq ) ) THEN 02158 DEALLOCATE ( pw_grid % gsq, STAT = stat ) 02159 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02160 END IF 02161 IF ( ASSOCIATED ( pw_grid % g_hat ) ) THEN 02162 DEALLOCATE ( pw_grid % g_hat, STAT = stat ) 02163 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02164 END IF 02165 IF ( ASSOCIATED ( pw_grid % grays ) ) THEN 02166 DEALLOCATE ( pw_grid % grays, STAT = stat ) 02167 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02168 END IF 02169 IF ( ASSOCIATED ( pw_grid % mapl % pos ) ) THEN 02170 DEALLOCATE ( pw_grid % mapl % pos, STAT = stat ) 02171 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02172 END IF 02173 IF ( ASSOCIATED ( pw_grid % mapm % pos ) ) THEN 02174 DEALLOCATE ( pw_grid % mapm % pos, STAT = stat ) 02175 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02176 END IF 02177 IF ( ASSOCIATED ( pw_grid % mapn % pos ) ) THEN 02178 DEALLOCATE ( pw_grid % mapn % pos, STAT = stat ) 02179 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02180 END IF 02181 IF ( ASSOCIATED ( pw_grid % mapl % neg ) ) THEN 02182 DEALLOCATE ( pw_grid % mapl % neg, STAT = stat ) 02183 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02184 END IF 02185 IF ( ASSOCIATED ( pw_grid % mapm % neg ) ) THEN 02186 DEALLOCATE ( pw_grid % mapm % neg, STAT = stat ) 02187 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02188 END IF 02189 IF ( ASSOCIATED ( pw_grid % mapn % neg ) ) THEN 02190 DEALLOCATE ( pw_grid % mapn % neg, STAT = stat ) 02191 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02192 END IF 02193 IF ( pw_grid % para % mode == PW_MODE_DISTRIBUTED ) THEN 02194 IF ( ASSOCIATED ( pw_grid % para % yzp ) ) THEN 02195 DEALLOCATE ( pw_grid % para % yzp, STAT = stat ) 02196 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02197 END IF 02198 IF ( ASSOCIATED ( pw_grid % para % yzq ) ) THEN 02199 DEALLOCATE ( pw_grid % para % yzq, STAT = stat ) 02200 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02201 END IF 02202 IF ( ASSOCIATED ( pw_grid % para % nyzray ) ) THEN 02203 DEALLOCATE ( pw_grid % para % nyzray, STAT = stat ) 02204 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02205 END IF 02206 IF ( ASSOCIATED ( pw_grid % para % bo ) ) THEN 02207 DEALLOCATE ( pw_grid % para % bo, STAT = stat ) 02208 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02209 END IF 02210 END IF 02211 ! also release groups 02212 CALL mp_comm_free ( pw_grid % para % group ) 02213 IF (PRODUCT(pw_grid % para % rs_dims) /= 0 ) & 02214 CALL mp_comm_free ( pw_grid % para % rs_group ) 02215 IF ( ASSOCIATED ( pw_grid % para % pos_of_x ) ) THEN 02216 DEALLOCATE ( pw_grid % para % pos_of_x, STAT = stat ) 02217 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02218 END IF 02219 02220 IF(ASSOCIATED(pw_grid)) THEN 02221 DEALLOCATE(pw_grid, stat=stat) 02222 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02223 END IF 02224 END IF 02225 END IF 02226 NULLIFY(pw_grid) 02227 END SUBROUTINE pw_grid_release 02228 02229 END MODULE pw_grids 02230
1.7.3