|
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 ! ***************************************************************************** 00013 MODULE se_ga_tools 00014 USE atomic_kind_types, ONLY: atomic_kind_type,& 00015 get_atomic_kind 00016 USE cell_types, ONLY: cell_type,& 00017 pbc 00018 USE cp_control_types, ONLY: dft_control_type,& 00019 semi_empirical_control_type 00020 USE cp_dbcsr_interface, ONLY: cp_dbcsr_iterator_blocks_left,& 00021 cp_dbcsr_iterator_next_block,& 00022 cp_dbcsr_iterator_start,& 00023 cp_dbcsr_iterator_stop 00024 USE cp_dbcsr_types, ONLY: cp_dbcsr_iterator,& 00025 cp_dbcsr_type 00026 USE cp_para_types, ONLY: cp_para_env_type 00027 USE distribution_1d_types, ONLY: distribution_1d_type 00028 USE ga_environment_types, ONLY: ga_environment_type 00029 USE kinds, ONLY: dp 00030 USE particle_types, ONLY: particle_type 00031 USE qs_environment_types, ONLY: get_qs_env,& 00032 qs_environment_type 00033 USE semi_empirical_types, ONLY: get_se_param,& 00034 semi_empirical_type 00035 USE timings, ONLY: timeset,& 00036 timestop 00037 #include "cp_common_uses.h" 00038 00039 IMPLICIT NONE 00040 PRIVATE 00041 00042 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'se_ga_tools' 00043 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. 00044 00045 00046 ! ***************************************************************************** 00047 00048 ! Public subroutine 00049 00050 PUBLIC :: se_deallocate_local_ga, se_ga_ks_accumulate, & 00051 se_allocate_local_ga, se_ga_initialize, & 00052 se_ga_allocate_local_info, se_ga_deallocate_local_info, & 00053 se_ga_get_nbin, se_ga_release, se_ga_diag_add, & 00054 se_ga_put_pmatrix, se_ga_pair_list_init 00055 00056 INTERFACE se_ga_allocate_local_info 00057 MODULE PROCEDURE se_ga_allocate_local_2, se_ga_allocate_local_1 00058 END INTERFACE 00059 00060 INTERFACE se_ga_deallocate_local_info 00061 MODULE PROCEDURE se_ga_deallocate_local_2, se_ga_deallocate_local_1 00062 END INTERFACE 00063 00064 INTERFACE se_ga_ks_accumulate 00065 MODULE PROCEDURE se_ga_ks_accumulate_1, se_ga_ks_accumulate_2 00066 END INTERFACE 00067 00068 00069 #ifdef _USE_GA 00070 #include "mafdecls.fh" 00071 #include "global.fh" 00072 #endif 00073 00074 CONTAINS 00075 ! ***************************************************************************** 00082 SUBROUTINE se_ga_release ( qs_env, error ) 00083 IMPLICIT NONE 00084 TYPE(qs_environment_type), POINTER :: qs_env 00085 TYPE(cp_error_type), INTENT(inout) :: error 00086 #ifdef _USE_GA 00087 INTEGER :: stat 00088 TYPE ( ga_environment_type ), POINTER :: ga_env 00089 LOGICAL :: status, failure 00090 CHARACTER(len=*), PARAMETER :: routineN = 'se_ga_release', 00091 routineP = moduleN//':'//routineN 00092 00093 CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure) 00094 NULLIFY ( ga_env ) 00095 00096 CALL get_qs_env(qs_env, ga_env=ga_env, error=error) 00097 00098 CPPrecondition(ASSOCIATED(ga_env),cp_failure_level,routineP,error,failure) 00099 00100 ! write timing 00101 ! WRITE ( *, * ) ga_nodeid (), ' TIME_PUT', ga_env % time_put, ga_env % iput 00102 ! WRITE ( *, * ) ga_nodeid (), ' TIME_GET', ga_env % time_get, ga_env % iget 00103 ! WRITE ( *, * ) ga_nodeid (), ' TIME_ACC', ga_env % time_acc, ga_env % iacc 00104 ! 00105 ! call ga_igop(1, ga_env % iput, 1, '+') 00106 ! call ga_igop(1, ga_env % iget, 1, '+') 00107 ! call ga_igop(1, ga_env % iacc, 1, '+') 00108 ! WRITE ( *,'(i4,a,i16)') ga_nodeid(),' Total puts: ', ga_env % iput 00109 ! WRITE ( *,'(i4,a,i16)') ga_nodeid(),' Total gets: ', ga_env % iget 00110 ! WRITE ( *,'(i4,a,i16)') ga_nodeid(),' Total accs: ', ga_env % iacc 00111 00112 ! call ga_igop(1, ga_env % ielm_get, 1, '+') 00113 ! call ga_igop(1, ga_env % ielm_put, 1, '+') 00114 ! call ga_igop(1, ga_env % ielm_acc, 1, '+') 00115 ! WRITE ( *,'(i4,a,i16)') ga_nodeid(),' Total elm puts: ', ga_env % ielm_put 00116 ! WRITE ( *,'(i4,a,i16)') ga_nodeid(),' Total elm gets: ', ga_env % ielm_get 00117 ! WRITE ( *,'(i4,a,i16)') ga_nodeid(),' Total elm accs: ', ga_env % ielm_acc 00118 00119 status=ga_destroy( ga_env % g_offset ) 00120 status=ga_destroy( ga_env % g_ks ) 00121 status=ga_destroy( ga_env % g_atoms ) 00122 status=ga_destroy( ga_env % g_d_info ) 00123 status=ga_destroy( ga_env % g_cellnum ) 00124 status=ga_destroy( ga_env % g_count ) 00125 status=ga_destroy( ga_env % g_cnt ) 00126 00127 #endif 00128 END SUBROUTINE se_ga_release 00129 ! ***************************************************************************** 00133 SUBROUTINE se_ga_pair_list_init ( qs_env, error ) 00134 IMPLICIT NONE 00135 TYPE(qs_environment_type), POINTER :: qs_env 00136 TYPE(cp_error_type), INTENT(inout) :: error 00137 #ifdef _USE_GA 00138 TYPE ( cp_para_env_type ), POINTER :: para_env 00139 TYPE ( ga_environment_type ), POINTER :: ga_env 00140 TYPE ( cell_type ), POINTER :: cell 00141 TYPE ( dft_control_type ), POINTER :: dft_control 00142 TYPE ( semi_empirical_control_type ), POINTER :: se_control 00143 REAL ( dp ) :: xbox, ybox, zbox 00144 REAL ( dp ) :: xcell, ycell, zcell, cx, cy, cz, rc, rcut2 00145 INTEGER :: kxmin,kymin,kzmin,kxmax,kymax,kzmax,kc,kx,ky,kz 00146 INTEGER :: nbins, nc, ix, iy, iz, itmp 00147 INTEGER :: jx, jy, jz, idx 00148 INTEGER :: ierr, ibeg, iend, me, nproc, i, j, icnt 00149 INTEGER :: g_cnt, one 00150 INTEGER :: lo(2), hi(2), dims(2),chunk(2), ndim, ld 00151 INTEGER :: pair(2) 00152 INTEGER :: pdx, pdy, pdz, number 00153 INTEGER :: kcmax, nix(2000), niy(2000), niz(2000) 00154 LOGICAL :: status, wrap 00155 00156 NULLIFY ( cell, ga_env, dft_control, se_control ) 00157 00158 CALL get_qs_env(qs_env, dft_control=dft_control, & 00159 para_env=para_env,ga_env=ga_env, & 00160 cell = cell, error=error ) 00161 00162 00163 se_control => dft_control % qs_control % se_control 00164 ga_env % rcoul = 2.0_dp * se_control % cutoff_cou 00165 ga_env % ncells = se_control % ga_ncells 00166 00167 WRITE ( *, * ) 'AAAA', ga_env % ncells, ga_env % rcoul 00168 00169 xbox=cell%hmat(1,1) 00170 ybox=cell%hmat(2,2) 00171 zbox=cell%hmat(3,3) 00172 00173 IF ( ga_env % ncells == 0 ) THEN 00174 ! use cut-off 00175 pdx=INT ( xbox/ga_env%rcoul ) 00176 IF ( ( xbox/REAL ( pdx, dp ) < ga_env % rcoul ) .AND. ( pdx > 1 ) ) THEN 00177 pdx = pdx - 1 00178 ELSE IF ( pdx == 0 ) THEN 00179 pdx = 1 00180 END IF 00181 pdy=INT ( ybox/ga_env%rcoul ) 00182 IF ( ( ybox/REAL ( pdy, dp ) < ga_env % rcoul ) .AND. ( pdy > 1 ) ) THEN 00183 pdy = pdy - 1 00184 ELSE IF ( pdy == 0 ) THEN 00185 pdy = 1 00186 END IF 00187 pdz=INT ( zbox/ga_env%rcoul ) 00188 IF ( ( zbox/REAL ( pdz, dp ) < ga_env % rcoul ) .AND. ( pdz > 1 ) ) THEN 00189 pdz = pdz - 1 00190 ELSE IF ( pdz == 0 ) THEN 00191 pdz = 1 00192 END IF 00193 ga_env % ncells = pdx*pdy*pdz 00194 WRITE ( *, * ) 'pds ', pdx,pdy,pdz 00195 ELSE 00196 ! use input driven number and compute pdx, pdy, pdz 00197 number = ga_env % ncells 00198 CALL grid_factor ( number, cell%hmat(1,1), cell%hmat(2,2), cell%hmat(3,3), pdx, pdy, pdz ) 00199 END IF 00200 ! Assign value for use later 00201 ga_env % pdx = pdx 00202 ga_env % pdy = pdy 00203 ga_env % pdz = pdz 00204 ! 00205 ! This finds the neighboring cells 00206 ! 00207 xcell = xbox / REAL ( pdx, dp ) 00208 kxmax = INT(ga_env % rcoul/xcell)+1 00209 ! 00210 ycell = ybox / REAL ( pdy, dp ) 00211 kymax = INT(ga_env % rcoul/ycell)+1 00212 ! 00213 zcell = zbox / REAL ( pdz, dp ) 00214 kzmax = INT(ga_env % rcoul/zcell)+1 00215 ! 00216 WRITE ( *, * ) 'ks ', kxmax, kymax, kzmax 00217 ! 00218 rcut2 = ga_env % rcoul**2 00219 ! 00220 kc = 1 00221 nix(kc) = 0 00222 niy(kc) = 0 00223 niz(kc) = 0 00224 ! 00225 kxmin = 0 00226 DO kx = kxmin, kxmax 00227 kymin = -kymax 00228 IF (kx==0) kymin = 0 00229 DO ky = kymin, kymax 00230 kzmin = -kzmax 00231 IF (kx==0.AND.ky==0) kzmin = 1 00232 DO kz = kzmin, kzmax 00233 ! write(6,'(a,6i3)') 'klims: ',kxmin,kxmax,kymin,kymax,kzmin,kzmax 00234 IF (kx/=0) THEN 00235 cx = xcell*REAL(iabs(kx)-1, dp) 00236 ELSE 00237 cx = 0.0_dp 00238 ENDIF 00239 IF (ky/=0) THEN 00240 cy = ycell*REAL(iabs(ky)-1, dp) 00241 ELSE 00242 cy = 0.0_dp 00243 ENDIF 00244 IF (kz/=0) THEN 00245 cz = zcell*REAL(iabs(kz)-1, dp) 00246 ELSE 00247 cz = 0.0_dp 00248 ENDIF 00249 rc = cx**2+cy**2+cz**2 00250 IF (rc<rcut2) THEN 00251 kc = kc + 1 00252 nix(kc) = kx 00253 niy(kc) = ky 00254 niz(kc) = kz 00255 ENDIF 00256 END DO 00257 END DO 00258 END DO 00259 kcmax = kc 00260 WRITE ( *, * ) 'kcmax ', kcmax 00261 ! We now have unique neighbors of each cell 00262 00263 ! 00264 ! Initialize system 00265 ! 00266 me = ga_nodeid() 00267 nproc = ga_nnodes() 00268 ! 00269 rcut2 = ga_env % rcoul**2 00270 xcell = xbox/REAL(pdx, dp) 00271 ycell = ybox/REAL(pdy, dp) 00272 zcell = zbox/REAL(pdz, dp) 00273 WRITE ( *, * ) 'boxs ', xbox, ybox, zbox 00274 WRITE ( *, * ) 'cells ', xcell, ycell, zcell 00275 ! 00276 ! Find total number of atom bins 00277 ! 00278 nbins = pdx*pdy*pdz 00279 WRITE ( *, * ) 'nbins ', nbins 00280 one = 1 00281 ! 00282 ! Find number of unique neighbors 00283 ! 00284 ! if (me.eq.0) write(6,'(a,i8)') 'Number of neighbors found: ',kcmax 00285 ! if (me.eq.0) then 00286 ! do i = 1, kcmax 00287 ! write(6,'(i4,a,3i4)') i,' nix, niy, niz: ',nix(i),niy(i),niz(i) 00288 ! end do 00289 ! endif 00290 ! 00291 ! Set up distributed pair list 00292 ! 00293 ibeg = INT(REAL(nbins,dp)*REAL(me,dp)/REAL(nproc,dp)) + 1 00294 IF (me.lt.nproc-1) THEN 00295 iend = INT(REAL(nbins,dp)*REAL(me+1,dp)/REAL(nproc,dp)) 00296 ELSE 00297 iend = nbins 00298 ENDIF 00299 ! 00300 ! Add all pairs to the distributed pair list 00301 ! 00302 icnt = 0 00303 DO i = ibeg, iend 00304 itmp = i - 1 00305 ix = MOD(itmp,pdx) 00306 itmp = (itmp-ix)/pdy 00307 iy = MOD(itmp,pdy) 00308 iz = (itmp-iy)/pdy 00309 ix = ix + 1 00310 iy = iy + 1 00311 iz = iz + 1 00312 DO nc = 1, kcmax 00313 wrap = .FALSE. 00314 jx = ix + nix(nc) 00315 IF (jx>pdx) THEN 00316 jx = MOD(jx,pdx) 00317 wrap = .TRUE. 00318 ENDIF 00319 IF (jx<1) THEN 00320 jx = pdx+jx 00321 wrap = .TRUE. 00322 ENDIF 00323 jy = iy + niy(nc) 00324 IF (jy>pdy) THEN 00325 jy = MOD(jy,pdy) 00326 wrap = .TRUE. 00327 ENDIF 00328 IF (jy<1) THEN 00329 jy = pdy+jy 00330 wrap = .TRUE. 00331 ENDIF 00332 jz = iz + niz(nc) 00333 IF (jz>pdz) THEN 00334 jz = MOD(jz,pdz) 00335 wrap = .TRUE. 00336 ENDIF 00337 IF (jz<1) THEN 00338 jz = pdz+jz 00339 wrap = .TRUE. 00340 ENDIF 00341 j = (jz-1)*pdx*pdy + (jy-1)*pdx + jx 00342 IF (wrap) THEN 00343 IF (iabs(ix-jx).ne.0) THEN 00344 cx = xcell*(iabs(ix-jx)-1) 00345 ELSE 00346 cx = 0.0_dp 00347 ENDIF 00348 IF (iabs(iy-jy)/=0) THEN 00349 cy = ycell*(iabs(iy-jy)-1) 00350 ELSE 00351 cy = 0.0_dp 00352 ENDIF 00353 IF (iabs(iz-jz)/=0) THEN 00354 cz = zcell*(iabs(iz-jz)-1) 00355 ELSE 00356 cz = 0.0_dp 00357 ENDIF 00358 rc = cx**2 + cy**2 + cz**2 00359 IF (rc<rcut2) CYCLE 00360 ENDIF 00361 icnt = icnt + 1 00362 END DO 00363 END DO 00364 ! 00365 ! Calculated all pair interactions on this processor 00366 ! 00367 ! Hold on to this 00368 ga_env % npairs = icnt 00369 WRITE ( *, * ) 'NPAIRS', icnt 00370 CALL ga_igop(1,ga_env%npairs,1,'+') 00371 ! if (me.eq.0) write(6,'(a,i8)') 'Total number of pairs found: ',npairs 00372 ! 00373 ! Create global counter and distributed pair list 00374 ! 00375 ndim = 1 00376 dims(1) = 1 00377 g_cnt = ga_create_handle() 00378 CALL ga_set_data(g_cnt, ndim, dims, MT_INT) 00379 status = ga_allocate(g_cnt) 00380 ! 00381 ndim = 2 00382 dims(1) = 2 00383 dims(2) = ga_env % npairs 00384 chunk(1) = 2 00385 chunk(2) = -1 00386 ! hold on g_pair_list 00387 ga_env % g_pair_list = ga_create_handle() 00388 CALL ga_set_data(ga_env % g_pair_list, ndim, dims, MT_INT) 00389 CALL ga_set_chunk(ga_env % g_pair_list, chunk) 00390 status = ga_allocate(ga_env % g_pair_list) 00391 ! 00392 CALL ga_zero(g_cnt) 00393 lo(1) = 1 00394 hi(1) = 2 00395 ld = 2 00396 DO i = ibeg, iend 00397 itmp = i - 1 00398 ix = MOD(itmp,pdx) 00399 itmp = (itmp-ix)/pdy 00400 iy = MOD(itmp,pdy) 00401 iz = (itmp-iy)/pdy 00402 ix = ix + 1 00403 iy = iy + 1 00404 iz = iz + 1 00405 pair(1) = i 00406 DO nc = 1, kcmax 00407 wrap = .FALSE. 00408 jx = ix + nix(nc) 00409 IF (jx>pdx) THEN 00410 jx = MOD(jx,pdx) 00411 wrap = .TRUE. 00412 ENDIF 00413 IF (jx<1) THEN 00414 jx = pdx+jx 00415 wrap = .TRUE. 00416 ENDIF 00417 jy = iy + niy(nc) 00418 IF (jy>pdy) THEN 00419 jy = MOD(jy,pdy) 00420 wrap = .TRUE. 00421 ENDIF 00422 IF (jy<1) THEN 00423 jy = pdy+jy 00424 wrap = .TRUE. 00425 ENDIF 00426 jz = iz + niz(nc) 00427 IF (jz>pdz) THEN 00428 jz = MOD(jz,pdz) 00429 wrap = .TRUE. 00430 ENDIF 00431 IF (jz<1) THEN 00432 jz = pdz+jz 00433 wrap = .TRUE. 00434 ENDIF 00435 j = (jz-1)*pdx*pdy + (jy-1)*pdx + jx 00436 IF (wrap) THEN 00437 IF (iabs(ix-jx)/=0) THEN 00438 cx = xcell*(iabs(ix-jx)-1) 00439 ELSE 00440 cx = 0.0_dp 00441 ENDIF 00442 IF (iabs(iy-jy)/=0) THEN 00443 cy = ycell*(iabs(iy-jy)-1) 00444 ELSE 00445 cy = 0.0_dp 00446 ENDIF 00447 IF (iabs(iz-jz)/=0) THEN 00448 cz = zcell*(iabs(iz-jz)-1) 00449 ELSE 00450 cz = 0.0_dp 00451 ENDIF 00452 rc = cx**2 + cy**2 + cz**2 00453 IF (rc<rcut2) CYCLE 00454 ENDIF 00455 idx = nga_read_inc(g_cnt, one, one) 00456 ! 00457 lo(2) = idx+1 00458 hi(2) = idx+1 00459 pair(2) = j 00460 CALL nga_put(ga_env % g_pair_list,lo,hi,pair,ld) 00461 END DO 00462 END DO 00463 CALL ga_sync 00464 ! 00465 ! Evaluate interactions between all pairs 00466 ! 00467 status = ga_destroy(g_cnt) 00468 #endif 00469 END SUBROUTINE se_ga_pair_list_init 00470 ! ***************************************************************************** 00474 SUBROUTINE se_ga_initialize ( qs_env, error ) 00475 IMPLICIT NONE 00476 TYPE(qs_environment_type), POINTER :: qs_env 00477 TYPE(cp_error_type), INTENT(inout) :: error 00478 #ifdef _USE_GA 00479 TYPE ( cp_para_env_type ), POINTER :: para_env 00480 TYPE ( ga_environment_type ), POINTER :: ga_env 00481 TYPE ( distribution_1d_type ), POINTER :: local_particles 00482 TYPE ( particle_type ), DIMENSION(:), 00483 POINTER :: particle_set 00484 TYPE ( cell_type ), POINTER :: cell 00485 TYPE ( semi_empirical_type ), POINTER :: se_kind 00486 TYPE ( atomic_kind_type ), DIMENSION(:), 00487 POINTER :: atomic_kind_set 00488 TYPE ( atomic_kind_type ), POINTER :: atomic_kind 00489 INTEGER :: nx, ny, nz, number, xbin, ybin, zbin, nbin 00490 INTEGER :: nparticle_kind, iparticle_kind, nparticle_local, iparticle_local 00491 INTEGER :: iparticle, k, j, icnt, natoms, natorb 00492 REAL ( dp ) :: dx, dy, dz, r_pbc ( 3 ) 00493 INTEGER :: stat, chunk ( 3 ), dims ( 3 ) 00494 INTEGER :: handle, lo ( 2 ), hi ( 2 ) 00495 INTEGER :: ioff, iloc 00496 INTEGER, PARAMETER :: one = 1 00497 INTEGER :: ii, jj 00498 INTEGER, SAVE :: max_rows 00499 LOGICAL, SAVE :: first_time=.TRUE. 00500 REAL ( dp ), POINTER :: block ( :, : ) 00501 REAL ( dp ), DIMENSION ( : ), POINTER :: tmp_info 00502 LOGICAL :: failure, status, check, defined, found 00503 CHARACTER(len=*), PARAMETER :: routineN = 'se_ga_intialize', 00504 routineP = moduleN//':'//routineN 00505 00506 CALL timeset(routineN,handle) 00507 CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure) 00508 NULLIFY ( para_env, cell, ga_env ) 00509 NULLIFY ( atomic_kind_set, particle_set, local_particles ) 00510 NULLIFY ( tmp_info ) 00511 00512 CALL get_qs_env(qs_env, para_env=para_env, & 00513 particle_set=particle_set, & 00514 atomic_kind_set=atomic_kind_set, & 00515 local_particles=local_particles, ga_env=ga_env, & 00516 cell = cell, error=error ) 00517 00518 00519 natoms = SIZE ( particle_set ) 00520 00521 00522 00523 CPPrecondition ( ASSOCIATED ( ga_env ), cp_failure_level,routineP,error,failure) 00524 00525 00526 number = ga_env % ncells 00527 nx = ga_env % pdx 00528 ny = ga_env % pdy 00529 nz = ga_env % pdz 00530 00531 ! CALL grid_factor ( number, cell%hmat(1,1), cell%hmat(2,2), cell%hmat(3,3), nx, ny, nz ) 00532 00533 dx=cell%hmat(1,1)/REAL ( nx, dp ) 00534 dy=cell%hmat(2,2)/REAL ( ny, dp ) 00535 dz=cell%hmat(3,3)/REAL ( nz, dp ) 00536 00537 ga_env % g_cnt=ga_create_handle( ) 00538 CALL ga_set_data ( ga_env % g_cnt, one, one, MT_F_INT ) 00539 status=ga_allocate ( ga_env % g_cnt ) 00540 CALL ga_zero ( ga_env % g_cnt ) 00541 00542 ga_env % g_cellnum=ga_create_handle( ) 00543 CALL ga_set_data ( ga_env % g_cellnum, one, number, MT_F_INT ) 00544 status=ga_allocate ( ga_env % g_cellnum ) 00545 CALL ga_zero ( ga_env % g_cellnum ) 00546 00547 ga_env % g_count=ga_create_handle ( ) 00548 CALL ga_set_data ( ga_env % g_count, one, number, MT_F_INT ) 00549 status=ga_allocate ( ga_env % g_count ) 00550 CALL ga_zero ( ga_env % g_count ) 00551 00552 ga_env % g_offset=ga_create_handle ( ) 00553 CALL ga_set_data ( ga_env % g_offset, one, number, MT_F_INT ) 00554 status=ga_allocate ( ga_env % g_offset ) 00555 CALL ga_zero ( ga_env % g_offset ) 00556 00557 IF ( para_env%mepos == 0 ) & 00558 CALL nga_put ( ga_env % g_count, one, one, one, one ) 00559 00560 ga_env % g_atoms=ga_create_handle ( ) 00561 CALL ga_set_data ( ga_env % g_atoms, 1, natoms, MT_F_INT ) 00562 status=ga_allocate ( ga_env % g_atoms ) 00563 CALL ga_zero ( ga_env % g_atoms ) 00564 00565 IF ( first_time ) max_rows = 0 00566 nparticle_kind = SIZE(atomic_kind_set) 00567 DO iparticle_kind=1,nparticle_kind 00568 atomic_kind => atomic_kind_set(iparticle_kind) 00569 CALL get_atomic_kind(atomic_kind=atomic_kind, se_parameter=se_kind ) 00570 CALL get_se_param ( se_kind, defined=defined, natorb=natorb ) 00571 IF (.NOT. defined .AND. natorb < 1) CYCLE 00572 nparticle_local = local_particles%n_el(iparticle_kind) 00573 DO iparticle_local=1,nparticle_local 00574 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 00575 IF ( first_time ) & 00576 max_rows = MAX ( natorb, max_rows ) 00577 r_pbc ( : ) = pbc ( particle_set ( iparticle ) % r ( : ), cell ) 00578 ! should be between 0 -> L 00579 CALL pbc_wrap ( r_pbc, cell%hmat(1,1), cell%hmat(2,2), cell%hmat (3,3) ) 00580 xbin = INT ( r_pbc ( 1 ) / dx ) 00581 ybin = INT ( r_pbc ( 2 ) / dy ) 00582 zbin = INT ( r_pbc ( 3 ) / dz ) 00583 nbin = zbin * nx * ny + ybin* nx + xbin + 1 00584 IF ( ( nbin > ga_env % ncells ) .OR. ( nbin < 1 ) ) THEN 00585 WRITE ( *, * ) "WARNING", iparticle, nbin, ga_env % ncells 00586 END IF 00587 CALL nga_acc ( ga_env % g_cellnum, nbin, nbin, one, one, one ) 00588 END DO 00589 END DO 00590 00591 IF ( first_time ) THEN 00592 ! CALL mp_max ( max_rows, para_env%group ) 00593 first_time = .FALSE. 00594 END IF 00595 !WRITE ( *, * ) 'MAXROWS', max_rows 00596 max_rows=4 00597 00598 ! set se_dims 00599 ga_env % se_dims ( 1 ) = max_rows 00600 ga_env % se_dims ( 2 ) = max_rows 00601 ga_env % se_dims ( 3 ) = natoms 00602 00603 ! compute offset array 00604 CALL ga_sync ( ) 00605 CALL ga_scan_add(ga_env % g_cellnum, ga_env % g_offset, ga_env % g_count, 1, number, 1) 00606 CALL ga_zero ( ga_env % g_count ) 00607 00608 dims(1) = ga_env % se_dims(1)*ga_env % se_dims(1) + 5 00609 dims(2) = ga_env % se_dims(3) 00610 chunk(1) = dims(1) 00611 chunk(2) = -1 00612 ALLOCATE(tmp_info(dims(1))) 00613 tmp_info = 0.0_dp 00614 00615 ga_env % g_d_info = ga_create_handle() 00616 CALL ga_set_data(ga_env % g_d_info,2,dims,MT_F_DBL) 00617 CALL ga_set_chunk(ga_env % g_d_info,chunk) 00618 status = ga_allocate(ga_env % g_d_info) 00619 00620 dims(1) = ga_env % se_dims(1) 00621 dims(2) = ga_env % se_dims(1) 00622 dims(3) = ga_env % se_dims(3) 00623 00624 chunk(1) = dims(1) 00625 chunk(2) = dims(1) 00626 chunk(3) = -1 00627 00628 ga_env % g_ks = ga_create_handle() 00629 CALL ga_set_data(ga_env % g_ks,3,dims,MT_F_DBL) 00630 CALL ga_set_chunk(ga_env % g_ks,chunk) 00631 status = ga_allocate(ga_env % g_ks) 00632 00633 00634 nparticle_kind = SIZE(atomic_kind_set) 00635 DO iparticle_kind=1,nparticle_kind 00636 atomic_kind => atomic_kind_set(iparticle_kind) 00637 CALL get_atomic_kind(atomic_kind=atomic_kind, se_parameter=se_kind ) 00638 CALL get_se_param ( se_kind, defined=defined, natorb=natorb ) 00639 ! WRITE ( *, * ) 'NATORB', natorb 00640 nparticle_local = local_particles%n_el(iparticle_kind) 00641 DO iparticle_local=1,nparticle_local 00642 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 00643 r_pbc ( : ) = pbc ( particle_set ( iparticle ) % r ( : ), cell ) 00644 ! should be between 0 -> L 00645 CALL pbc_wrap ( r_pbc, cell%hmat(1,1), cell%hmat(2,2), cell%hmat (3,3) ) 00646 xbin = INT ( r_pbc ( 1 ) / dx ) 00647 ybin = INT ( r_pbc ( 2 ) / dy ) 00648 zbin = INT ( r_pbc ( 3 ) / dz ) 00649 nbin = zbin * nx * ny + ybin* nx + xbin + 1 00650 iloc = nga_read_inc(ga_env % g_count, nbin, 1) 00651 CALL nga_get(ga_env % g_offset, nbin, nbin, ioff, 1) 00652 CALL nga_put(ga_env % g_atoms, iparticle, iparticle, ioff+iloc+1,1) 00653 00654 tmp_info(1) = REAL ( iparticle_kind, dp ) 00655 tmp_info(2) = REAL ( iparticle, dp ) 00656 tmp_info(3) = particle_set ( iparticle ) % r ( 1 ) 00657 tmp_info(4) = particle_set ( iparticle ) % r ( 2 ) 00658 tmp_info(5) = particle_set ( iparticle ) % r ( 3 ) 00659 ! WRITE ( *, * ) 'TMP ', iparticle, SNGL ( tmp_info ( 3:5 )) 00660 icnt = 0 00661 DO k = 1, natorb 00662 DO j = 1, natorb 00663 icnt = icnt + 1 00664 tmp_info(5+icnt) = 0.0_dp 00665 END DO 00666 END DO 00667 lo(1) = 1 00668 lo(2) = ioff+iloc+1 00669 hi(1) = ga_env%se_dims ( 1 )**2 + 5 00670 hi(2) = ioff+iloc+1 00671 CALL nga_put(ga_env % g_d_info, lo, hi, tmp_info, dims(1)) 00672 END DO 00673 END DO 00674 DEALLOCATE(tmp_info) 00675 CALL ga_sync () 00676 CALL timestop(handle) 00677 #endif 00678 END SUBROUTINE se_ga_initialize 00679 ! ***************************************************************************** 00686 SUBROUTINE se_allocate_local_ga ( buff, max_rows, error ) 00687 REAL ( dp ), POINTER :: buff ( :, : ) 00688 INTEGER, INTENT ( IN ) :: max_rows 00689 TYPE(cp_error_type), INTENT(inout) :: error 00690 #ifdef _USE_GA 00691 INTEGER :: stat 00692 LOGICAL :: failure 00693 CHARACTER(len=*), PARAMETER :: routineN = 'se_allocate_local_ga', 00694 routineP = moduleN//':'//routineN 00695 ALLOCATE ( buff ( max_rows, max_rows ), STAT=stat ) 00696 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00697 #endif 00698 END SUBROUTINE se_allocate_local_ga 00699 ! ***************************************************************************** 00706 SUBROUTINE se_deallocate_local_ga ( buff, error ) 00707 REAL ( dp ), POINTER :: buff ( :, : ) 00708 TYPE(cp_error_type), INTENT(inout) :: error 00709 #ifdef _USE_GA 00710 INTEGER :: stat 00711 LOGICAL :: failure, status 00712 CHARACTER(len=*), PARAMETER :: routineN = 'se_deallocate_local_ga', 00713 routineP = moduleN//':'//routineN 00714 DEALLOCATE(buff,STAT=stat) 00715 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00716 NULLIFY ( buff ) 00717 #endif 00718 END SUBROUTINE se_deallocate_local_ga 00719 ! ***************************************************************************** 00720 SUBROUTINE se_ga_put_pmatrix ( matrix, ga_env, error ) 00721 IMPLICIT NONE 00722 TYPE ( cp_dbcsr_type ), POINTER :: matrix 00723 TYPE ( ga_environment_type ), POINTER :: ga_env 00724 TYPE(cp_error_type), INTENT(inout) :: error 00725 #ifdef _USE_GA 00726 TYPE ( cp_dbcsr_iterator ) :: iter 00727 LOGICAL :: failure 00728 INTEGER :: handle, stat 00729 INTEGER :: iblock_row, iblock_col, rows_block 00730 INTEGER :: lo ( 2 ), hi ( 2 ), ld 00731 INTEGER :: blk, idx 00732 INTEGER, PARAMETER :: one = 1 00733 REAL(kind=dp), DIMENSION(:, :), POINTER :: block 00734 CHARACTER(len=*), PARAMETER :: routineN = 'se_ga_put_pmatrix', 00735 routineP = moduleN//':'//routineN 00736 00737 CALL timeset(routineN,handle) 00738 CALL ga_sync 00739 CALL cp_dbcsr_iterator_start (iter, matrix ) 00740 DO WHILE (cp_dbcsr_iterator_blocks_left(iter)) 00741 CALL cp_dbcsr_iterator_next_block (iter, iblock_row,iblock_col, block ,blk) 00742 IF ( iblock_row == iblock_col) THEN 00743 rows_block=SIZE(block,1) 00744 CALL nga_get ( ga_env % g_atoms, iblock_row, iblock_row, idx, one ) 00745 lo(1) = 6 00746 lo(2) = idx 00747 hi(1) = 5+rows_block**2 00748 hi(2) = idx 00749 ld = 1 00750 CALL nga_put(ga_env % g_d_info, lo, hi, block, ld) 00751 END IF 00752 END DO 00753 CALL cp_dbcsr_iterator_stop ( iter ) 00754 CALL timestop(handle) 00755 #endif 00756 END SUBROUTINE se_ga_put_pmatrix 00757 ! ***************************************************************************** 00758 SUBROUTINE se_ga_diag_add ( matrix, ga_env, error ) 00759 IMPLICIT NONE 00760 TYPE ( cp_dbcsr_type ), POINTER :: matrix 00761 TYPE ( ga_environment_type ), POINTER :: ga_env 00762 TYPE(cp_error_type), INTENT(inout) :: error 00763 #ifdef _USE_GA 00764 TYPE ( cp_dbcsr_iterator ) :: iter 00765 LOGICAL :: failure 00766 INTEGER :: handle, stat 00767 INTEGER :: iblock_row, iblock_col, rows_block 00768 INTEGER :: lo ( 3 ), hi ( 3 ), ld ( 2 ) 00769 INTEGER :: blk, idx 00770 INTEGER, PARAMETER :: one = 1 00771 REAL(kind=dp), DIMENSION(:, :), POINTER :: dum, block 00772 CHARACTER(len=*), PARAMETER :: routineN = 'se_ga_diag_add', 00773 routineP = moduleN//':'//routineN 00774 00775 CALL timeset(routineN,handle) 00776 CALL ga_sync 00777 CALL cp_dbcsr_iterator_start (iter, matrix ) 00778 DO WHILE (cp_dbcsr_iterator_blocks_left(iter)) 00779 CALL cp_dbcsr_iterator_next_block (iter, iblock_row,iblock_col, block ,blk) 00780 IF ( iblock_row == iblock_col) THEN 00781 rows_block=SIZE(block,1) 00782 CALL nga_get ( ga_env % g_atoms, iblock_row, iblock_row, idx, one ) 00783 lo(1) = 1 00784 lo(2) = 1 00785 lo(3) = idx 00786 hi(1) = rows_block 00787 hi(2) = rows_block 00788 hi(3) = idx 00789 ld(1) = rows_block 00790 ld(2) = rows_block 00791 ALLOCATE ( dum ( rows_block, rows_block ),STAT=stat ) 00792 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00793 CALL nga_get(ga_env % g_ks, lo, hi, dum, ld) 00794 block = block + dum 00795 DEALLOCATE ( dum , STAT=stat) 00796 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00797 END IF 00798 END DO 00799 CALL cp_dbcsr_iterator_stop ( iter ) 00800 CALL timestop(handle) 00801 #endif 00802 END SUBROUTINE se_ga_diag_add 00803 ! ***************************************************************************** 00804 SUBROUTINE se_ga_ks_accumulate_2 ( ga_env, ks_i, ks_j, ioff, joff, isize, jsize ) 00805 TYPE ( ga_environment_type ), POINTER :: ga_env 00806 REAL ( dp ), DIMENSION ( :, :, : ), POINTER :: ks_i, ks_j 00807 INTEGER, INTENT ( IN ) :: ioff, joff, isize, jsize 00808 #ifdef _USE_GA 00809 INTEGER :: ilo ( 3 ), jlo ( 3 ), ihi ( 3 ), jhi ( 3 ) 00810 REAL ( dp ), PARAMETER :: one=1.0_dp 00811 ilo(1) = 1 00812 ilo(2) = 1 00813 ilo(3) = ioff + 1 00814 jlo(1) = 1 00815 jlo(2) = 1 00816 jlo(3) = joff + 1 00817 ihi(1) = ga_env % se_dims(1) 00818 ihi(2) = ga_env % se_dims(2) 00819 ihi(3) = ioff + isize 00820 jhi(1) = ga_env % se_dims(1) 00821 jhi(2) = ga_env % se_dims(2) 00822 jhi(3) = joff + jsize 00823 CALL nga_acc ( ga_env % g_ks, ilo, ihi, ks_i, ga_env % se_dims ( 1 ), one ) 00824 CALL nga_acc ( ga_env % g_ks, jlo, jhi, ks_j, ga_env % se_dims ( 1 ), one ) 00825 #endif 00826 END SUBROUTINE se_ga_ks_accumulate_2 00827 ! ***************************************************************************** 00828 SUBROUTINE se_ga_ks_accumulate_1 ( ga_env, ks_i, ioff, isize ) 00829 TYPE ( ga_environment_type ), POINTER :: ga_env 00830 REAL ( dp ), DIMENSION ( :, :, : ), POINTER :: ks_i 00831 INTEGER, INTENT ( IN ) :: ioff, isize 00832 #ifdef _USE_GA 00833 INTEGER :: ilo ( 3 ), ihi ( 3 ) 00834 REAL ( dp ), PARAMETER :: one=1.0_dp 00835 ilo(1) = 1 00836 ilo(2) = 1 00837 ilo(3) = ioff + 1 00838 ihi(1) = ga_env % se_dims(1) 00839 ihi(2) = ga_env % se_dims(2) 00840 ihi(3) = ioff + isize 00841 CALL nga_acc ( ga_env % g_ks, ilo, ihi, ks_i, ga_env % se_dims ( 1 ), one ) 00842 #endif 00843 END SUBROUTINE se_ga_ks_accumulate_1 00844 ! ***************************************************************************** 00845 SUBROUTINE se_ga_deallocate_local_2 ( dens_i_info, dens_j_info, ks_i, ks_j, error ) 00846 REAL ( dp ), POINTER :: dens_i_info ( :, : ), dens_j_info ( :, : ) 00847 REAL ( dp ), POINTER :: ks_i ( :, :, : ), ks_j ( :, :, : ) 00848 TYPE(cp_error_type), INTENT(inout) :: error 00849 #ifdef _USE_GA 00850 INTEGER :: stat 00851 CHARACTER(len=*), PARAMETER :: routineN = 'se_ga_deallocate_local_2', 00852 routineP = moduleN//':'//routineN 00853 LOGICAL :: failure 00854 00855 IF ( ASSOCIATED ( dens_i_info ) ) DEALLOCATE ( dens_i_info, STAT=stat ) 00856 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00857 IF ( ASSOCIATED ( dens_j_info ) ) DEALLOCATE ( dens_j_info, STAT=stat ) 00858 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00859 IF ( ASSOCIATED ( ks_i ) ) DEALLOCATE ( ks_i, STAT=stat ) 00860 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00861 IF ( ASSOCIATED ( ks_j ) ) DEALLOCATE ( ks_j, STAT=stat ) 00862 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00863 NULLIFY ( dens_i_info ) 00864 NULLIFY ( dens_j_info ) 00865 NULLIFY ( ks_i ) 00866 NULLIFY ( ks_j ) 00867 #endif 00868 END SUBROUTINE se_ga_deallocate_local_2 00869 ! ***************************************************************************** 00870 SUBROUTINE se_ga_deallocate_local_1 ( dens_i_info, ks_i, error ) 00871 REAL ( dp ), POINTER, OPTIONAL :: dens_i_info ( :, : ) 00872 REAL ( dp ), POINTER, OPTIONAL:: ks_i ( :, :, : ) 00873 TYPE(cp_error_type), INTENT(inout) :: error 00874 #ifdef _USE_GA 00875 INTEGER :: stat 00876 CHARACTER(len=*), PARAMETER :: routineN = 'se_ga_deallocate_local_1', 00877 routineP = moduleN//':'//routineN 00878 LOGICAL :: failure 00879 00880 IF ( PRESENT ( dens_i_info ) ) THEN 00881 IF ( ASSOCIATED ( dens_i_info ) ) DEALLOCATE ( dens_i_info, STAT=stat ) 00882 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00883 NULLIFY ( dens_i_info ) 00884 END IF 00885 IF ( PRESENT ( ks_i ) ) THEN 00886 IF ( ASSOCIATED ( ks_i ) ) DEALLOCATE ( ks_i, STAT=stat ) 00887 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00888 NULLIFY ( ks_i ) 00889 END IF 00890 #endif 00891 END SUBROUTINE se_ga_deallocate_local_1 00892 ! ***************************************************************************** 00893 SUBROUTINE se_ga_allocate_local_2 ( ga_env, dens_i_info, dens_j_info, ks_i, ks_j, & 00894 ibin, jbin, isize, jsize, ioff, joff, error ) 00895 TYPE ( ga_environment_type ), POINTER :: ga_env 00896 INTEGER, INTENT ( IN ) :: ibin, jbin 00897 INTEGER, INTENT ( OUT ) :: isize, jsize, ioff, joff 00898 REAL ( dp ), POINTER :: dens_i_info ( :, : ), dens_j_info ( :, : ) 00899 REAL ( dp ), POINTER :: ks_i ( :, :, : ), ks_j ( :, :, : ) 00900 TYPE(cp_error_type), INTENT(inout) :: error 00901 #ifdef _USE_GA 00902 CHARACTER(len=*), PARAMETER :: routineN = 'se_ga_allocate_local_2', 00903 routineP = moduleN//':'//routineN 00904 LOGICAL :: failure 00905 INTEGER, PARAMETER :: one=1 00906 INTEGER :: ilo ( 2 ), jlo ( 2 ), ihi ( 2 ), jhi ( 2 ), stat 00907 CALL nga_get(ga_env%g_offset,ibin,ibin,ioff,one) 00908 CALL nga_get(ga_env%g_cellnum,ibin,ibin,isize,one) 00909 CALL nga_get(ga_env%g_offset,jbin,jbin,joff,one) 00910 CALL nga_get(ga_env%g_cellnum,jbin,jbin,jsize,one) 00911 IF ( ( isize == 0 ) .OR. ( jsize == 0 ) ) RETURN 00912 ilo(1) = 1 00913 ilo(2) = ioff + 1 00914 jlo(1) = 1 00915 jlo(2) = joff + 1 00916 ihi(1) = ga_env % se_dims(1)*ga_env % se_dims(1) + 5 00917 ihi(2) = ioff + isize 00918 jhi(1) = ga_env % se_dims(1)*ga_env % se_dims(1) + 5 00919 jhi(2) = joff + jsize 00920 ALLOCATE ( dens_i_info ( ihi ( 1 ), isize ), STAT=stat ) 00921 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00922 ALLOCATE ( dens_j_info ( jhi ( 1 ), jsize ), STAT=stat ) 00923 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00924 ALLOCATE ( ks_i ( ga_env % se_dims ( 1 ), ga_env % se_dims ( 1 ), isize ), STAT=stat ) 00925 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00926 ALLOCATE ( ks_j ( ga_env % se_dims ( 1 ), ga_env % se_dims ( 1 ), jsize ), STAT=stat ) 00927 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00928 ks_i = 0.0_dp 00929 ks_j = 0.0_dp 00930 CALL nga_get(ga_env % g_d_info,ilo,ihi,dens_i_info,ihi(1)) 00931 CALL nga_get(ga_env % g_d_info,jlo,jhi,dens_j_info,jhi(1)) 00932 #else 00933 isize = 1 00934 jsize = 1 00935 ioff = 0 00936 joff = 0 00937 #endif 00938 END SUBROUTINE se_ga_allocate_local_2 00939 ! ***************************************************************************** 00940 SUBROUTINE se_ga_allocate_local_1 ( ga_env, dens_i_info, ks_i,& 00941 ibin, isize, ioff, error ) 00942 TYPE ( ga_environment_type ), POINTER :: ga_env 00943 INTEGER, INTENT ( IN ) :: ibin 00944 INTEGER, INTENT ( OUT ):: isize 00945 INTEGER, INTENT ( OUT ), OPTIONAL :: ioff 00946 REAL ( dp ), POINTER :: dens_i_info ( :, : ) 00947 REAL ( dp ), POINTER, OPTIONAL :: ks_i ( :, :, : ) 00948 TYPE(cp_error_type), INTENT(inout) :: error 00949 #ifdef _USE_GA 00950 CHARACTER(len=*), PARAMETER :: routineN = 'se_ga_allocate_local_1', 00951 routineP = moduleN//':'//routineN 00952 LOGICAL :: failure 00953 INTEGER :: myioff 00954 INTEGER, PARAMETER :: one=1 00955 INTEGER :: ilo ( 2 ), ihi ( 2 ), stat 00956 CALL nga_get(ga_env%g_offset,ibin,ibin,myioff,one) 00957 CALL nga_get(ga_env%g_cellnum,ibin,ibin,isize,one) 00958 IF ( isize == 0 ) RETURN 00959 IF ( PRESENT ( ioff ) ) ioff = myioff 00960 ilo(1) = 1 00961 ilo(2) = myioff + 1 00962 ihi(1) = ga_env % se_dims(1) * ga_env % se_dims(1) + 5 00963 ihi(2) = myioff + isize 00964 ALLOCATE ( dens_i_info ( ihi ( 1 ), isize ), STAT=stat ) 00965 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00966 CALL nga_get(ga_env % g_d_info,ilo,ihi,dens_i_info,ihi(1)) 00967 IF ( PRESENT ( ks_i ) ) THEN 00968 ALLOCATE ( ks_i ( ga_env % se_dims ( 1 ), ga_env % se_dims ( 1 ), isize ), STAT=stat ) 00969 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00970 ks_i = 0.0_dp 00971 END IF 00972 #else 00973 isize = 1 00974 IF ( PRESENT ( ioff ) ) ioff = 0 00975 #endif 00976 END SUBROUTINE se_ga_allocate_local_1 00977 00978 ! ***************************************************************************** 00979 SUBROUTINE se_ga_get_nbin ( ga_handle, nbin, lzero ) 00980 INTEGER, INTENT ( IN ) :: ga_handle 00981 INTEGER, INTENT ( OUT ) :: nbin 00982 LOGICAL, INTENT ( IN ), OPTIONAL :: lzero 00983 #ifdef _USE_GA 00984 LOGICAL :: my_lzero 00985 INTEGER, PARAMETER :: one=1 00986 my_lzero = .FALSE. 00987 IF ( PRESENT ( lzero ) ) my_lzero = lzero 00988 IF ( my_lzero ) CALL ga_zero ( ga_handle ) 00989 nbin = nga_read_inc ( ga_handle, one, one ) 00990 #else 00991 nbin = 1 00992 #endif 00993 END SUBROUTINE se_ga_get_nbin 00994 00995 ! ***************************************************************************** 00996 SUBROUTINE grid_factor(p,xdim,ydim,zdim,idx,idy,idz) 00997 INTEGER, INTENT(IN) :: p 00998 REAL(dp), INTENT(IN) :: xdim, ydim, zdim 00999 INTEGER, INTENT(OUT) :: idx, idy, idz 01000 01001 INTEGER, PARAMETER :: MAX_FACTOR = 10000 01002 01003 INTEGER :: fac(MAX_FACTOR), i, ifac, ip, 01004 ix, iy, iz, j, pmax, 01005 prime(MAX_FACTOR) 01006 01007 ! 01008 01009 i = 1 01010 ! 01011 ! factor p completely 01012 ! first, find all prime numbers, besides 1, less than or equal to 01013 ! the square root of p 01014 ! 01015 ip = INT(SQRT(DBLE(p)))+1 01016 pmax = 0 01017 DO i = 2, ip 01018 DO j = 1, pmax 01019 IF (MOD(i,prime(j))==0) go to 100 01020 END DO 01021 pmax = pmax + 1 01022 IF (pmax>MAX_FACTOR) WRITE(6,*) 'Overflow in grid_factor' 01023 prime(pmax) = i 01024 100 CONTINUE 01025 END DO 01026 ! 01027 ! find all prime factors of p 01028 ! 01029 ip = p 01030 ifac = 0 01031 DO i = 1, pmax 01032 200 IF (MOD(ip,prime(i))==0) THEN 01033 ifac = ifac + 1 01034 fac(ifac) = prime(i) 01035 ip = ip/prime(i) 01036 go to 200 01037 ENDIF 01038 END DO 01039 ! 01040 ! determine three factors of p of approximately the 01041 ! same size 01042 ! 01043 idx = 1 01044 idy = 1 01045 idz = 1 01046 DO i = ifac, 1, -1 01047 ix = xdim / idx 01048 iy = ydim / idy 01049 iz = zdim / idz 01050 IF (ix>=iy.and.ix>=iz.and.ix>1) THEN 01051 idx = fac(i)*idx 01052 ELSEIF (iy>=ix.and.iy>=iz.and.iy>1) THEN 01053 idy = fac(i)*idy 01054 ELSEIF (iz>=ix.and.iz>=iy.and.iz>1) THEN 01055 idz = fac(i)*idz 01056 ELSE 01057 WRITE(6,*) 'Too many processors in grid factoring routine' 01058 ENDIF 01059 END DO 01060 END SUBROUTINE grid_factor 01061 ! ***************************************************************************** 01062 SUBROUTINE pbc_wrap ( r, x, y, z ) 01063 REAL(kind=dp), DIMENSION(3), 01064 INTENT(INOUT) :: r 01065 REAL(kind=dp), INTENT(IN) :: x, y, z 01066 01067 DO WHILE ( r ( 1 ) > x ) 01068 r ( 1 ) = r ( 1 ) - x 01069 END DO 01070 DO WHILE ( r ( 2 ) > y ) 01071 r ( 2 ) = r ( 2 ) - y 01072 END DO 01073 DO WHILE ( r ( 3 ) > z ) 01074 r ( 3 ) = r ( 3 ) - z 01075 END DO 01076 DO WHILE ( r ( 1 ) < 0.0_dp ) 01077 r ( 1 ) = r ( 1 ) + x 01078 END DO 01079 DO WHILE ( r ( 2 ) < 0.0_dp ) 01080 r ( 2 ) = r ( 2 ) + y 01081 END DO 01082 DO WHILE ( r ( 3 ) < 0.0_dp ) 01083 r ( 3 ) = r ( 3 ) + z 01084 END DO 01085 END SUBROUTINE pbc_wrap 01086 ! ***************************************************************************** 01087 END MODULE se_ga_tools
1.7.3