|
CP2K 2.4 (Revision 12889)
|
00001 00002 !-----------------------------------------------------------------------------! 00003 ! CP2K: A general program to perform molecular dynamics simulations ! 00004 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00005 !-----------------------------------------------------------------------------! 00006 00007 ! ***************************************************************************** 00023 MODULE pw_methods 00024 USE f77_blas 00025 USE fft_tools, ONLY: BWFFT,& 00026 FWFFT,& 00027 fft3d 00028 USE kahan_sum, ONLY: accurate_sum 00029 USE kinds, ONLY: dp 00030 USE machine, ONLY: m_loc_c,& 00031 m_loc_r 00032 USE message_passing, ONLY: mp_sum 00033 USE pw_grid_types, ONLY: HALFSPACE,& 00034 PW_MODE_DISTRIBUTED,& 00035 PW_MODE_LOCAL,& 00036 pw_grid_type 00037 USE pw_methods_cuda, ONLY: cuda_pw_fft_wrap_c1dr3d,& 00038 cuda_pw_fft_wrap_r3dc1d 00039 USE pw_types, ONLY: & 00040 COMPLEXDATA1D, COMPLEXDATA3D, NOSPACE, REALDATA1D, REALDATA3D, & 00041 REALSPACE, RECIPROCALSPACE, SQUARE, SQUAREROOT, pw_type 00042 USE termination, ONLY: stop_program 00043 USE timings, ONLY: timeset,& 00044 timestop 00045 #include "cp_common_uses.h" 00046 00047 IMPLICIT NONE 00048 00049 PRIVATE 00050 00051 PUBLIC :: pw_zero, pw_structure_factor, pw_smoothing 00052 PUBLIC :: pw_copy, pw_axpy, pw_transfer, pw_scale 00053 PUBLIC :: pw_derive, pw_dr2, pw_fft_wrap, pw_write 00054 PUBLIC :: pw_compare_debug 00055 PUBLIC :: pw_integral_aa, pw_integral_ab, pw_integral_a2b 00056 PUBLIC :: pw_integral_aa_fast, pw_integral_ab_fast 00057 PUBLIC :: pw_dr2_gg, pw_integrate_function 00058 00059 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_methods' 00060 LOGICAL, PARAMETER, PRIVATE :: debug_this_module=.FALSE. 00061 00062 INTERFACE pw_gather 00063 MODULE PROCEDURE pw_gather_s, pw_gather_p 00064 END INTERFACE 00065 00066 INTERFACE pw_scatter 00067 MODULE PROCEDURE pw_scatter_s, pw_scatter_p 00068 END INTERFACE 00069 00070 INTERFACE pw_fft_wrap 00071 MODULE PROCEDURE fft_wrap_pw1, fft_wrap_pw1pw2 00072 END INTERFACE 00073 00074 CONTAINS 00075 00076 ! ***************************************************************************** 00082 SUBROUTINE pw_zero ( pw, error ) 00083 00084 TYPE(pw_type), INTENT(INOUT) :: pw 00085 TYPE(cp_error_type), INTENT(inout) :: error 00086 00087 CHARACTER(len=*), PARAMETER :: routineN = 'pw_zero', 00088 routineP = moduleN//':'//routineN 00089 00090 INTEGER :: handle, ns 00091 LOGICAL :: failure 00092 REAL(KIND=dp) :: zr 00093 00094 failure = .FALSE. 00095 CALL timeset(routineN,handle) 00096 CPPrecondition(pw%ref_count>0,cp_failure_level,routineP,error,failure) 00097 IF ( pw%in_use == REALDATA1D ) THEN 00098 ns = SIZE(pw%cr) 00099 pw%cr(:) = 0._dp 00100 ELSE IF ( pw%in_use == COMPLEXDATA1D ) THEN 00101 ns = SIZE(pw%cc) 00102 pw%cc(:) = CMPLX(0._dp,0._dp,KIND=dp) 00103 ELSE IF ( pw%in_use == REALDATA3D ) THEN 00104 ns = SIZE(pw%cr3d) 00105 pw%cr3d(:,:,:) = 0._dp 00106 ELSE IF ( pw%in_use == COMPLEXDATA3D ) THEN 00107 ns = SIZE(pw%cc3d) 00108 pw%cc3d(:,:,:) = CMPLX(0._dp,0._dp,KIND=dp) 00109 ELSE 00110 CALL stop_program(routineN,moduleN,__LINE__,"No possible data field!") 00111 END IF 00112 00113 zr = REAL ( ns,KIND=dp) * 1.e-6_dp 00114 CALL timestop(handle) 00115 00116 END SUBROUTINE pw_zero 00117 00118 ! ***************************************************************************** 00129 SUBROUTINE pw_copy ( pw1, pw2, error) 00130 TYPE(pw_type), INTENT(IN) :: pw1 00131 TYPE(pw_type), INTENT(INOUT) :: pw2 00132 TYPE(cp_error_type), INTENT(inout) :: error 00133 00134 CHARACTER(len=*), PARAMETER :: routineN = 'pw_copy', 00135 routineP = moduleN//':'//routineN 00136 00137 INTEGER :: handle, i, j, ng, ng1, ng2, ns 00138 LOGICAL :: failure 00139 REAL(KIND=dp) :: zc 00140 00141 failure = .FALSE. 00142 CALL timeset(routineN,handle) 00143 CPPrecondition(pw1%ref_count>0,cp_failure_level,routineP,error,failure) 00144 CPPrecondition(pw2%ref_count>0,cp_failure_level,routineP,error,failure) 00145 IF ( pw1%pw_grid %id_nr /= pw2%pw_grid %id_nr ) THEN 00146 00147 IF ( pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical ) THEN 00148 00149 IF ( pw_compatible (pw1%pw_grid, pw2%pw_grid, error=error) ) THEN 00150 00151 IF ( pw1%in_use == COMPLEXDATA1D .AND. & 00152 pw2%in_use == COMPLEXDATA1D .AND. & 00153 pw1%in_space == RECIPROCALSPACE ) THEN 00154 ng1 = SIZE ( pw1%cc ) 00155 ng2 = SIZE ( pw2%cc ) 00156 ng = MIN ( ng1, ng2 ) 00157 !$omp parallel do private(i) 00158 DO i = 1, ng 00159 pw2%cc(i) = pw1%cc(i) 00160 END DO 00161 IF ( ng2 > ng ) THEN 00162 !$omp parallel do private(i) 00163 DO i = ng+1, ng2 00164 pw2%cc(i) = CMPLX ( 0.0_dp, 0.0_dp,KIND=dp) 00165 END DO 00166 END IF 00167 ns = 2 * MAX ( ng1, ng2 ) 00168 ELSE 00169 CALL stop_program(routineN,moduleN,__LINE__,"No suitable data field") 00170 END IF 00171 00172 ELSE 00173 WRITE ( *, "(A,I5,T30,A,L1,T60,A,I5)" ) & 00174 " grid 1 :",pw1%pw_grid %id_nr, & 00175 " spherical :",pw1%pw_grid%spherical, & 00176 " reference :",pw1%pw_grid%reference 00177 WRITE ( *, "(A,I5,T30,A,L1,T60,A,I5)" ) & 00178 " grid 2 :",pw2%pw_grid %id_nr, & 00179 " spherical :",pw2%pw_grid%spherical, & 00180 " reference :",pw2%pw_grid%reference 00181 CALL stop_program(routineN,moduleN,__LINE__,"Incompatible grids") 00182 END IF 00183 00184 ELSE IF ( .NOT. ( pw1%pw_grid%spherical .OR. & 00185 pw2%pw_grid%spherical ) ) THEN 00186 00187 ng1 = SIZE ( pw1%cc ) 00188 ng2 = SIZE ( pw2%cc ) 00189 ns = 2 * MAX ( ng1, ng2 ) 00190 00191 IF ( pw1%in_use == COMPLEXDATA1D .AND. & 00192 pw2%in_use == COMPLEXDATA1D .AND. & 00193 pw1%in_space == RECIPROCALSPACE ) THEN 00194 00195 IF ( ( pw1%pw_grid %id_nr == pw2%pw_grid%reference ) ) THEN 00196 IF( ng1 >= ng2 ) THEN 00197 !$omp parallel do private(i,j) 00198 DO i = 1, ng2 00199 j = pw2%pw_grid%gidx ( i ) 00200 pw2%cc ( i ) = pw1%cc ( j ) 00201 END DO 00202 ELSE 00203 CALL pw_zero ( pw2, error=error) 00204 !$omp parallel do private(i,j) 00205 DO i = 1, ng1 00206 j = pw2%pw_grid%gidx ( i ) 00207 pw2%cc ( j ) = pw1%cc ( i ) 00208 END DO 00209 END IF 00210 ELSE IF ( ( pw2%pw_grid %id_nr == pw1%pw_grid%reference ) ) THEN 00211 IF( ng1 >= ng2 ) THEN 00212 !$omp parallel do private(i,j) 00213 DO i = 1, ng2 00214 j = pw1%pw_grid%gidx ( i ) 00215 pw2%cc ( i ) = pw1%cc ( j ) 00216 END DO 00217 ELSE 00218 CALL pw_zero ( pw2, error=error) 00219 !$omp parallel do private(i,j) 00220 DO i = 1, ng1 00221 j = pw1%pw_grid%gidx ( i ) 00222 pw2%cc ( j ) = pw1%cc ( i ) 00223 END DO 00224 END IF 00225 ELSE 00226 WRITE ( *, "(A,I5,T30,A,L1,T60,A,I5)" ) & 00227 " grid 1 :",pw1%pw_grid %id_nr, & 00228 " spherical :",pw1%pw_grid%spherical, & 00229 " reference :",pw1%pw_grid%reference 00230 WRITE ( *, "(A,I5,T30,A,L1,T60,A,I5)" ) & 00231 " grid 2 :",pw2%pw_grid %id_nr, & 00232 " spherical :",pw2%pw_grid%spherical, & 00233 " reference :",pw2%pw_grid%reference 00234 CALL stop_program(routineN,moduleN,__LINE__,"Incompatible grids") 00235 END IF 00236 00237 ELSE 00238 CALL stop_program(routineN,moduleN,__LINE__,"No suitable data field") 00239 END IF 00240 00241 pw2%in_space = RECIPROCALSPACE 00242 00243 ELSE 00244 WRITE ( *, "(A,I5,T30,A,L1,T60,A,I5)" ) & 00245 " grid 1 :",pw1%pw_grid %id_nr, & 00246 " spherical :",pw1%pw_grid%spherical, & 00247 " reference :",pw1%pw_grid%reference 00248 WRITE ( *, "(A,I5,T30,A,L1,T60,A,I5)" ) & 00249 " grid 2 :",pw2%pw_grid %id_nr, & 00250 " spherical :",pw2%pw_grid%spherical, & 00251 " reference :",pw2%pw_grid%reference 00252 CALL stop_program(routineN,moduleN,__LINE__,"Incompatible grids") 00253 END IF 00254 00255 ELSE 00256 00257 IF ( pw1%in_use == REALDATA1D .AND. pw2%in_use == REALDATA1D ) THEN 00258 ns = SIZE ( pw1%cr ) 00259 !$omp parallel do private(i) 00260 DO i = 1, ns 00261 pw2%cr(i) = pw1%cr(i) 00262 END DO 00263 ELSE IF ( pw1%in_use == COMPLEXDATA1D .AND. & 00264 pw2%in_use == COMPLEXDATA1D ) THEN 00265 ns = SIZE ( pw1%cc ) 00266 !$omp parallel do private(i) 00267 DO i = 1, ns 00268 pw2%cc(i) = pw1%cc(i) 00269 END DO 00270 ELSE IF ( pw1%in_use == REALDATA3D .AND. pw2%in_use == REALDATA3D ) THEN 00271 ns = SIZE ( pw1%cr3d ) 00272 pw2%cr3d(:,:,:) = pw1%cr3d(:,:,:) 00273 ELSE IF ( pw1%in_use == COMPLEXDATA3D .AND. & 00274 pw2%in_use == COMPLEXDATA3D ) THEN 00275 ns = SIZE ( pw1%cc3d ) 00276 pw2%cc3d(:,:,:) = pw1%cc3d(:,:,:) 00277 ELSE 00278 CALL stop_program(routineN,moduleN,__LINE__,"No suitable data field") 00279 END IF 00280 00281 END IF 00282 00283 pw2%in_space = pw1%in_space 00284 zc = REAL ( ns,KIND=dp) * 1.e-6_dp 00285 CALL timestop(handle) 00286 00287 END SUBROUTINE pw_copy 00288 00289 ! ***************************************************************************** 00294 SUBROUTINE pw_scale( pw, a , error) 00295 00296 TYPE(pw_type), INTENT(INOUT) :: pw 00297 REAL(KIND=dp), INTENT(IN) :: a 00298 TYPE(cp_error_type), INTENT(inout) :: error 00299 00300 CHARACTER(len=*), PARAMETER :: routineN = 'pw_scale', 00301 routineP = moduleN//':'//routineN 00302 00303 INTEGER :: handle, ns 00304 LOGICAL :: failure 00305 REAL(KIND=dp) :: flop 00306 00307 failure = .FALSE. 00308 CALL timeset(routineN,handle) 00309 CPPrecondition(pw%ref_count>0,cp_failure_level,routineP,error,failure) 00310 00311 SELECT CASE ( pw%in_use ) 00312 CASE ( REALDATA1D ) 00313 ns = SIZE(pw%cr) 00314 pw%cr(:) = a*pw%cr(:) 00315 CASE ( COMPLEXDATA1D ) 00316 ns = 2*SIZE(pw%cc) 00317 pw%cc(:) = a*pw%cc(:) 00318 CASE ( REALDATA3D) 00319 ns = SIZE(pw%cr3d) 00320 pw%cr3d(:,:,:) = a*pw%cr3d(:,:,:) 00321 CASE (COMPLEXDATA3D ) 00322 ns = 2*SIZE(pw%cc3d) 00323 pw%cc3d(:,:,:) = a*pw%cc3d(:,:,:) 00324 CASE DEFAULT 00325 CALL stop_program(routineN,moduleN,__LINE__, "No suitable data field" ) 00326 END SELECT 00327 00328 flop = REAL ( ns, KIND=dp) * 1.e-6_dp 00329 CALL timestop(handle) 00330 00331 END SUBROUTINE pw_scale 00332 00333 ! ***************************************************************************** 00342 SUBROUTINE pw_derive ( pw, n, error ) 00343 00344 TYPE(pw_type), INTENT(INOUT) :: pw 00345 INTEGER, DIMENSION(3), INTENT(IN) :: n 00346 TYPE(cp_error_type), INTENT(inout) :: error 00347 00348 CHARACTER(len=*), PARAMETER :: routineN = 'pw_derive', 00349 routineP = moduleN//':'//routineN 00350 00351 COMPLEX(KIND=dp) :: im 00352 INTEGER :: cnt, handle, i, m 00353 LOGICAL :: failure 00354 REAL(KIND=dp) :: flop 00355 00356 failure = .FALSE. 00357 CALL timeset(routineN,handle) 00358 CPPrecondition(pw%ref_count>0,cp_failure_level,routineP,error,failure) 00359 CPPrecondition(ALL(n>=0),cp_failure_level,routineP,error,failure) 00360 00361 m = SUM ( n ) 00362 im = CMPLX ( 0.0_dp, 1.0_dp,KIND=dp) ** m 00363 00364 flop = 0.0_dp 00365 00366 IF ( pw%in_space == RECIPROCALSPACE .AND. & 00367 pw%in_use == COMPLEXDATA1D ) THEN 00368 00369 cnt = SIZE ( pw%cc ) 00370 00371 IF ( n ( 1 ) == 1 ) THEN 00372 !$omp parallel do private (i) 00373 DO i = 1, cnt 00374 pw%cc ( i ) = pw%cc ( i ) * pw%pw_grid%g ( 1, i ) 00375 END DO 00376 flop = flop + 6 * cnt 00377 ELSE IF ( n ( 1 ) > 1 ) THEN 00378 !$omp parallel do private (i) 00379 DO i = 1, cnt 00380 pw%cc ( i ) = pw%cc ( i ) * ( pw%pw_grid%g ( 1, i ) ** n ( 1 ) ) 00381 END DO 00382 flop = flop + 7 * cnt 00383 END IF 00384 IF ( n ( 2 ) == 1 ) THEN 00385 !$omp parallel do private (i) 00386 DO i = 1, cnt 00387 pw%cc ( i ) = pw%cc ( i ) * pw%pw_grid%g ( 2, i ) 00388 END DO 00389 flop = flop + 6 * cnt 00390 ELSE IF ( n ( 2 ) > 1 ) THEN 00391 !$omp parallel do private (i) 00392 DO i = 1, cnt 00393 pw%cc ( i ) = pw%cc ( i ) * ( pw%pw_grid%g ( 2, i ) ** n ( 2 ) ) 00394 END DO 00395 flop = flop + 7 * cnt 00396 END IF 00397 IF ( n ( 3 ) == 1 ) THEN 00398 !$omp parallel do private (i) 00399 DO i = 1, cnt 00400 pw%cc ( i ) = pw%cc ( i ) * pw%pw_grid%g ( 3, i ) 00401 END DO 00402 flop = flop + 6 * cnt 00403 ELSE IF ( n ( 3 ) > 1 ) THEN 00404 !$omp parallel do private (i) 00405 DO i = 1, cnt 00406 pw%cc ( i ) = pw%cc ( i ) * ( pw%pw_grid%g ( 3, i ) ** n ( 3 ) ) 00407 END DO 00408 flop = flop + 7 * cnt 00409 END IF 00410 00411 ! im can take the values 1, -1, i, -i 00412 ! skip this if im == 1 00413 IF ( ABS ( REAL ( im,KIND=dp) - 1.0_dp ) > 1.e-10 ) THEN 00414 !$omp parallel do private (i) 00415 DO i = 1, cnt 00416 pw%cc ( i ) = im * pw%cc ( i ) 00417 END DO 00418 flop = flop + 6 * cnt 00419 END IF 00420 00421 ELSE 00422 00423 CALL stop_program(routineN,moduleN,__LINE__,"No suitable data field") 00424 00425 END IF 00426 00427 flop = flop * 1.e-6_dp 00428 CALL timestop(handle) 00429 00430 END SUBROUTINE pw_derive 00431 00432 ! ***************************************************************************** 00440 SUBROUTINE pw_dr2 ( pw, pwdr2, i, j, error) 00441 00442 TYPE(pw_type), INTENT(IN) :: pw 00443 TYPE(pw_type), INTENT(INOUT) :: pwdr2 00444 INTEGER, INTENT(IN) :: i, j 00445 TYPE(cp_error_type), INTENT(inout) :: error 00446 00447 CHARACTER(len=*), PARAMETER :: routineN = 'pw_dr2', 00448 routineP = moduleN//':'//routineN 00449 00450 INTEGER :: cnt, handle, ig 00451 LOGICAL :: failure 00452 REAL(KIND=dp) :: flop, gg, o3 00453 00454 failure = .FALSE. 00455 CALL timeset(routineN,handle) 00456 CPPrecondition(pw%ref_count>0,cp_failure_level,routineP,error,failure) 00457 00458 flop = 0.0_dp 00459 o3 = 1._dp/3._dp 00460 00461 IF ( pw%in_space == RECIPROCALSPACE .AND. & 00462 pw%in_use == COMPLEXDATA1D ) THEN 00463 00464 cnt = SIZE ( pw%cc ) 00465 00466 IF ( i == j ) THEN 00467 !$omp parallel do private (ig,gg) default(none) shared(i,o3,pw,pwdr2,cnt) 00468 DO ig = 1, cnt 00469 gg = pw%pw_grid%g(i,ig)**2 - o3*pw%pw_grid%gsq(ig) 00470 pwdr2%cc(ig) = gg * pw%cc(ig) 00471 END DO 00472 flop = flop + 5 * cnt 00473 ELSE 00474 !$omp parallel do private (ig) 00475 DO ig = 1, cnt 00476 pwdr2%cc(ig) = pw%cc(ig) * (pw%pw_grid%g(i,ig)*pw%pw_grid%g(j,ig)) 00477 END DO 00478 flop = flop + 4 * cnt 00479 END IF 00480 00481 ELSE 00482 00483 CALL stop_program(routineN,moduleN,__LINE__,"No suitable data field") 00484 00485 END IF 00486 00487 flop = flop * 1.e-6_dp 00488 CALL timestop(handle) 00489 00490 END SUBROUTINE pw_dr2 00491 00492 ! ***************************************************************************** 00502 SUBROUTINE pw_dr2_gg ( pw, pwdr2_gg, i, j, error) 00503 00504 TYPE(pw_type), INTENT(IN) :: pw 00505 TYPE(pw_type), INTENT(INOUT) :: pwdr2_gg 00506 INTEGER, INTENT(IN) :: i, j 00507 TYPE(cp_error_type), INTENT(inout) :: error 00508 00509 INTEGER :: cnt, handle, ig 00510 LOGICAL :: failure 00511 REAL(KIND=dp) :: flop, gg, o3 00512 CHARACTER(len=*), PARAMETER :: routineN = 'pw_dr2_gg', 00513 routineP = moduleN//':'//routineN 00514 00515 failure = .FALSE. 00516 CALL timeset(routineN,handle) 00517 CPPrecondition(pw%ref_count>0,cp_failure_level,routineP,error,failure) 00518 00519 flop = 0.0_dp 00520 o3 = 1._dp/3._dp 00521 00522 IF ( pw%in_space == RECIPROCALSPACE .AND. & 00523 pw%in_use == COMPLEXDATA1D ) THEN 00524 00525 cnt = SIZE ( pw%cc ) 00526 00527 IF ( i == j ) THEN 00528 !$omp parallel do private (ig) 00529 DO ig = pw%pw_grid%first_gne0, cnt 00530 gg = pw%pw_grid%g(i,ig)**2 - o3*pw%pw_grid%gsq(ig) 00531 pwdr2_gg%cc(ig) = gg * pw%cc(ig) / pw%pw_grid%gsq(ig) 00532 END DO 00533 flop = flop + 6 * cnt 00534 ELSE 00535 !$omp parallel do private (ig) 00536 DO ig = pw%pw_grid%first_gne0, cnt 00537 pwdr2_gg%cc(ig) = pw%cc(ig) * (pw%pw_grid%g(i,ig)*pw%pw_grid%g(j,ig)) & 00538 / pw%pw_grid%gsq(ig) 00539 END DO 00540 flop = flop + 5 * cnt 00541 END IF 00542 00543 IF ( pw%pw_grid%have_g0 ) pwdr2_gg%cc ( 1 ) = 0.0_dp 00544 00545 ELSE 00546 00547 CALL stop_program(routineN,moduleN,__LINE__,"No suitable data field") 00548 00549 END IF 00550 00551 flop = flop * 1.e-6_dp 00552 CALL timestop(handle) 00553 00554 END SUBROUTINE pw_dr2_gg 00555 00556 ! ***************************************************************************** 00565 SUBROUTINE pw_smoothing ( pw, ecut, sigma, error) 00566 00567 TYPE(pw_type), INTENT(INOUT) :: pw 00568 REAL(KIND=dp), INTENT(IN) :: ecut, sigma 00569 TYPE(cp_error_type), INTENT(inout) :: error 00570 00571 CHARACTER(len=*), PARAMETER :: routineN = 'pw_smoothing', 00572 routineP = moduleN//':'//routineN 00573 00574 INTEGER :: cnt, handle, ig 00575 LOGICAL :: failure 00576 REAL(KIND=dp) :: arg, f, flop 00577 00578 failure = .FALSE. 00579 CALL timeset(routineN,handle) 00580 CPPrecondition(pw%ref_count>0,cp_failure_level,routineP,error,failure) 00581 00582 flop = 0.0_dp 00583 00584 IF ( pw%in_space == RECIPROCALSPACE .AND. & 00585 pw%in_use == COMPLEXDATA1D ) THEN 00586 00587 cnt = SIZE ( pw%cc ) 00588 00589 !$omp parallel do private (ig,arg,f) default(none) shared(sigma,ecut,pw,cnt) 00590 DO ig = 1, cnt 00591 arg = (ecut - pw%pw_grid%gsq(ig))/sigma 00592 f = EXP(arg)/(1+EXP(arg)) 00593 pw%cc(ig) = f * pw%cc(ig) 00594 END DO 00595 flop = flop + 6 * cnt 00596 00597 ELSE 00598 00599 CALL stop_program(routineN,moduleN,__LINE__,"No suitable data field") 00600 00601 END IF 00602 00603 flop = flop * 1.e-6_dp 00604 CALL timestop(handle) 00605 00606 END SUBROUTINE pw_smoothing 00607 00608 ! ***************************************************************************** 00616 SUBROUTINE pw_transfer ( pw1, pw2, debug, error) 00617 00618 TYPE(pw_type), INTENT(IN), TARGET :: pw1 00619 TYPE(pw_type), INTENT(INOUT), TARGET :: pw2 00620 LOGICAL, INTENT(IN), OPTIONAL :: debug 00621 TYPE(cp_error_type), INTENT(inout) :: error 00622 00623 CHARACTER(len=*), PARAMETER :: routineN = 'pw_transfer', 00624 routineP = moduleN//':'//routineN 00625 00626 LOGICAL :: failure 00627 00628 failure = .FALSE. 00629 CPPrecondition(pw1%ref_count>0,cp_failure_level,routineP,error,failure) 00630 CPPrecondition(pw2%ref_count>0,cp_failure_level,routineP,error,failure) 00631 00632 IF ( pw1%in_space == REALSPACE .AND. pw2%in_space == REALSPACE ) THEN 00633 00634 ! simple copy should do 00635 CALL pw_copy ( pw1, pw2, error=error) 00636 00637 ELSEIF ( pw1%in_space == RECIPROCALSPACE .AND. & 00638 pw2%in_space == RECIPROCALSPACE ) THEN 00639 00640 IF ( pw1%in_use == pw2%in_use ) THEN 00641 00642 ! simple copy should do 00643 CALL pw_copy ( pw1, pw2, error=error) 00644 00645 ELSE 00646 00647 ! we have to gather/scatter the data 00648 IF ( pw1%in_use == COMPLEXDATA1D ) THEN 00649 CALL pw_scatter ( pw1, pw2%cc3d, error=error) 00650 ELSEIF ( pw2%in_use == COMPLEXDATA1D ) THEN 00651 CALL pw_gather ( pw2, pw1%cc3d, error=error) 00652 ELSE 00653 CALL stop_program(routineN,moduleN,__LINE__, "Do not know what to do" ) 00654 END IF 00655 00656 END IF 00657 00658 ELSE 00659 00660 ! FFT needed, all further tests done in pw_fft_wrap 00661 CALL pw_fft_wrap ( pw1, pw2, debug, error=error) 00662 00663 END IF 00664 00665 END SUBROUTINE pw_transfer 00666 00667 ! ***************************************************************************** 00678 SUBROUTINE pw_axpy ( pw1, pw2, alpha, error) 00679 00680 TYPE(pw_type), INTENT(IN) :: pw1 00681 TYPE(pw_type), INTENT(INOUT) :: pw2 00682 REAL(KIND=dp), INTENT(in), OPTIONAL :: alpha 00683 TYPE(cp_error_type), INTENT(inout) :: error 00684 00685 CHARACTER(len=*), PARAMETER :: routineN = 'pw_axpy', 00686 routineP = moduleN//':'//routineN 00687 00688 INTEGER :: handle, i, j, ng, ng1, ng2 00689 LOGICAL :: failure 00690 REAL(KIND=dp) :: flop, my_alpha 00691 00692 failure = .FALSE. 00693 CALL timeset(routineN,handle) 00694 CPPrecondition(pw1%ref_count>0,cp_failure_level,routineP,error,failure) 00695 CPPrecondition(pw2%ref_count>0,cp_failure_level,routineP,error,failure) 00696 00697 my_alpha=1.0_dp 00698 IF (PRESENT(alpha)) my_alpha=alpha 00699 00700 IF ( pw1%pw_grid %id_nr == pw2%pw_grid %id_nr ) THEN 00701 00702 IF ( pw1%in_use == REALDATA1D .AND. pw2%in_use == REALDATA1D ) THEN 00703 IF (my_alpha==1.0_dp) THEN 00704 !$omp parallel do private(i) 00705 DO i = 1, SIZE ( pw2%cr ) 00706 pw2%cr(i) = pw2%cr(i) + pw1%cr(i) 00707 END DO 00708 flop = REAL ( SIZE ( pw2%cr ),KIND=dp) 00709 ELSE 00710 !$omp parallel do private(i) 00711 DO i = 1, SIZE ( pw2%cr ) 00712 pw2%cr(i) = pw2%cr(i) + my_alpha*pw1%cr(i) 00713 END DO 00714 flop = REAL ( 2*SIZE ( pw2%cr ),KIND=dp) 00715 END IF 00716 ELSE IF ( pw1%in_use == COMPLEXDATA1D .AND. & 00717 pw2%in_use == COMPLEXDATA1D ) THEN 00718 IF (my_alpha==1.0_dp) THEN 00719 !$omp parallel do private(i) 00720 DO i = 1, SIZE ( pw2%cc ) 00721 pw2%cc(i) = pw2%cc(i) + pw1%cc(i) 00722 END DO 00723 flop = REAL ( 2 * SIZE ( pw2%cc ),KIND=dp) 00724 ELSE 00725 !$omp parallel do private(i) 00726 DO i = 1, SIZE ( pw2%cc ) 00727 pw2%cc(i) = pw2%cc(i) + my_alpha*pw1%cc(i) 00728 END DO 00729 flop = REAL ( 4 * SIZE ( pw2%cc ),KIND=dp) 00730 END IF 00731 ELSE IF ( pw1%in_use == REALDATA3D .AND. pw2%in_use == REALDATA3D ) THEN 00732 IF (my_alpha==1.0_dp) THEN 00733 pw2%cr3d = pw2%cr3d + pw1%cr3d 00734 flop = REAL ( SIZE ( pw2%cr3d ),KIND=dp) 00735 ELSE 00736 pw2%cr3d = pw2%cr3d + my_alpha * pw1%cr3d 00737 flop = REAL ( 2*SIZE ( pw2%cr3d ),KIND=dp) 00738 END IF 00739 ELSE IF ( pw1%in_use == COMPLEXDATA3D .AND. & 00740 pw2%in_use == COMPLEXDATA3D ) THEN 00741 IF (my_alpha==1.0_dp) THEN 00742 pw2%cc3d = pw2%cc3d + pw1%cc3d 00743 flop = REAL ( 2 * SIZE ( pw2%cc3d ),KIND=dp) 00744 ELSE 00745 pw2%cc3d = pw2%cc3d + my_alpha*pw1%cc3d 00746 flop = REAL ( 4 * SIZE ( pw2%cc3d ),KIND=dp) 00747 END IF 00748 ELSE 00749 CALL stop_program(routineN,moduleN,__LINE__,"No suitable data field") 00750 END IF 00751 00752 ELSE IF ( pw_compatible ( pw1%pw_grid, pw2%pw_grid, error=error) ) THEN 00753 00754 IF ( pw1%in_use == COMPLEXDATA1D .AND. & 00755 pw2%in_use == COMPLEXDATA1D .AND. & 00756 pw1%in_space == RECIPROCALSPACE .AND. & 00757 pw1%in_space == RECIPROCALSPACE ) THEN 00758 00759 ng1 = SIZE ( pw1%cc ) 00760 ng2 = SIZE ( pw2%cc ) 00761 ng = MIN ( ng1, ng2 ) 00762 flop = REAL ( 2 * ng,KIND=dp) 00763 00764 IF ( pw1%pw_grid%spherical ) THEN 00765 00766 IF (my_alpha==1.0_dp) THEN 00767 !$omp parallel do private(i) 00768 DO i = 1, ng 00769 pw2%cc (i) = pw2%cc (i) + pw1%cc (i) 00770 END DO 00771 ELSE 00772 !$omp parallel do private(i) 00773 DO i = 1, ng 00774 pw2%cc (i) = pw2%cc (i) + my_alpha*pw1%cc (i) 00775 END DO 00776 END IF 00777 00778 ELSEIF ( ( pw1%pw_grid %id_nr == pw2%pw_grid%reference ) ) THEN 00779 00780 IF (my_alpha==1.0_dp) THEN 00781 IF( ng1 >= ng2 ) THEN 00782 !$omp parallel do private(i,j) 00783 DO i = 1, ng2 00784 j = pw2%pw_grid%gidx ( i ) 00785 pw2%cc ( i ) = pw2%cc ( i ) + pw1%cc ( j ) 00786 END DO 00787 ELSE 00788 !$omp parallel do private(i,j) 00789 DO i = 1, ng1 00790 j = pw2%pw_grid%gidx ( i ) 00791 pw2%cc ( j ) = pw2%cc ( j ) + pw1%cc ( i ) 00792 END DO 00793 END IF 00794 ELSE 00795 IF( ng1 >= ng2 ) THEN 00796 !$omp parallel do private(i,j) 00797 DO i = 1, ng2 00798 j = pw2%pw_grid%gidx ( i ) 00799 pw2%cc ( i ) = pw2%cc ( i ) + my_alpha*pw1%cc ( j ) 00800 END DO 00801 ELSE 00802 !$omp parallel do private(i,j) 00803 DO i = 1, ng1 00804 j = pw2%pw_grid%gidx ( i ) 00805 pw2%cc ( j ) = pw2%cc ( j ) + my_alpha*pw1%cc ( i ) 00806 END DO 00807 END IF 00808 END IF 00809 ELSEIF ( ( pw2%pw_grid %id_nr == pw1%pw_grid%reference ) ) THEN 00810 00811 IF (my_alpha==1.0_dp) THEN 00812 IF( ng1 >= ng2 ) THEN 00813 !$omp parallel do private(i,j) 00814 DO i = 1, ng2 00815 j = pw1%pw_grid%gidx ( i ) 00816 pw2%cc ( i ) = pw2%cc ( i ) + pw1%cc ( j ) 00817 END DO 00818 ELSE 00819 !$omp parallel do private(i,j) 00820 DO i = 1, ng1 00821 j = pw1%pw_grid%gidx ( i ) 00822 pw2%cc ( j ) = pw2%cc ( j ) + pw1%cc ( i ) 00823 END DO 00824 END IF 00825 ELSE 00826 IF( ng1 >= ng2 ) THEN 00827 !$omp parallel do private(i,j) 00828 DO i = 1, ng2 00829 j = pw1%pw_grid%gidx ( i ) 00830 pw2%cc ( i ) = pw2%cc ( i ) + my_alpha* pw1%cc ( j ) 00831 END DO 00832 ELSE 00833 !$omp parallel do private(i,j) 00834 DO i = 1, ng1 00835 j = pw1%pw_grid%gidx ( i ) 00836 pw2%cc ( j ) = pw2%cc ( j ) + my_alpha*pw1%cc ( i ) 00837 END DO 00838 END IF 00839 END IF 00840 ELSE 00841 00842 WRITE ( *, "(A,I5,T30,A,L1,T60,A,I5)" ) & 00843 " grid 1 :",pw1%pw_grid %id_nr, & 00844 " sperical :",pw1%pw_grid%spherical, & 00845 " reference :",pw1%pw_grid%reference 00846 WRITE ( *, "(A,I5,T30,A,L1,T60,A,I5)" ) & 00847 " grid 2 :",pw2%pw_grid %id_nr, & 00848 " sperical :",pw2%pw_grid%spherical, & 00849 " reference :",pw2%pw_grid%reference 00850 CALL stop_program(routineN,moduleN,__LINE__,"Grids not compatible") 00851 00852 END IF 00853 00854 ELSE 00855 CALL stop_program(routineN,moduleN,__LINE__,"No suitable data field") 00856 END IF 00857 00858 ELSE 00859 00860 WRITE ( *, "(A,I5,T30,A,L1,T60,A,I5)" ) & 00861 " grid 1 :",pw1%pw_grid %id_nr, & 00862 " sperical :",pw1%pw_grid%spherical, & 00863 " reference :",pw1%pw_grid%reference 00864 WRITE ( *, "(A,I5,T30,A,L1,T60,A,I5)" ) & 00865 " grid 2 :",pw2%pw_grid %id_nr, & 00866 " sperical :",pw2%pw_grid%spherical, & 00867 " reference :",pw2%pw_grid%reference 00868 CALL stop_program(routineN,moduleN,__LINE__,"Grids not compatible") 00869 00870 END IF 00871 00872 flop = flop * 1.e-6_dp 00873 CALL timestop(handle) 00874 00875 END SUBROUTINE pw_axpy 00876 00877 ! ***************************************************************************** 00883 SUBROUTINE pw_gather_s ( pw, c, scale, error) 00884 00885 TYPE(pw_type), INTENT(INOUT) :: pw 00886 COMPLEX(KIND=dp), DIMENSION(:, :, :), 00887 INTENT(IN) :: c 00888 REAL(KIND=dp), INTENT(IN), OPTIONAL :: scale 00889 TYPE(cp_error_type), INTENT(inout) :: error 00890 00891 CHARACTER(len=*), PARAMETER :: routineN = 'pw_gather_s', 00892 routineP = moduleN//':'//routineN 00893 00894 INTEGER :: gpt, handle, l, m, n, ngpts 00895 INTEGER, DIMENSION(:), POINTER :: mapl, mapm, mapn 00896 INTEGER, DIMENSION(:, :), POINTER :: ghat 00897 LOGICAL :: failure 00898 REAL(KIND=dp) :: cpy 00899 00900 failure = .FALSE. 00901 CALL timeset(routineN,handle) 00902 CPPrecondition(pw%ref_count>0,cp_failure_level,routineP,error,failure) 00903 00904 IF ( pw%in_use /= COMPLEXDATA1D ) THEN 00905 CALL stop_program(routineN,moduleN,__LINE__,"Data field has to be COMPLEXDATA1D") 00906 ENDIF 00907 00908 ! after the gather we are in g-space 00909 pw%in_space = RECIPROCALSPACE 00910 00911 mapl => pw%pw_grid%mapl%pos 00912 mapm => pw%pw_grid%mapm%pos 00913 mapn => pw%pw_grid%mapn%pos 00914 00915 ngpts = SIZE ( pw%pw_grid%gsq ) 00916 00917 ghat => pw%pw_grid%g_hat 00918 00919 IF ( PRESENT ( scale ) ) THEN 00920 00921 !$omp parallel do private(gpt,l,m,n) 00922 DO gpt = 1, ngpts 00923 00924 l = mapl ( ghat ( 1, gpt ) ) + 1 00925 m = mapm ( ghat ( 2, gpt ) ) + 1 00926 n = mapn ( ghat ( 3, gpt ) ) + 1 00927 pw%cc ( gpt ) = scale * c ( l, m, n ) 00928 00929 END DO 00930 00931 ELSE 00932 00933 !$omp parallel do private(gpt,l,m,n) 00934 DO gpt = 1, ngpts 00935 00936 l = mapl ( ghat ( 1, gpt ) ) + 1 00937 m = mapm ( ghat ( 2, gpt ) ) + 1 00938 n = mapn ( ghat ( 3, gpt ) ) + 1 00939 pw%cc ( gpt ) = c ( l, m, n ) 00940 00941 END DO 00942 00943 END IF 00944 00945 cpy = REAL ( ngpts,KIND=dp) * 1.e-6_dp 00946 CALL timestop(handle) 00947 00948 END SUBROUTINE pw_gather_s 00949 00950 ! ***************************************************************************** 00951 SUBROUTINE pw_gather_p ( pw, c, scale, error) 00952 00953 TYPE(pw_type), INTENT(INOUT), TARGET :: pw 00954 COMPLEX(KIND=dp), DIMENSION(:, :), 00955 INTENT(IN) :: c 00956 REAL(KIND=dp), INTENT(IN), OPTIONAL :: scale 00957 TYPE(cp_error_type), INTENT(inout) :: error 00958 00959 CHARACTER(len=*), PARAMETER :: routineN = 'pw_gather_p', 00960 routineP = moduleN//':'//routineN 00961 00962 INTEGER :: gpt, handle, l, m, mn, n, 00963 ngpts 00964 INTEGER, DIMENSION(:), POINTER :: mapl, mapm, mapn 00965 INTEGER, DIMENSION(:, :), POINTER :: ghat, yzq 00966 LOGICAL :: failure 00967 REAL(KIND=dp) :: cpy 00968 00969 failure = .FALSE. 00970 CALL timeset(routineN,handle) 00971 CPPrecondition(pw%ref_count>0,cp_failure_level,routineP,error,failure) 00972 00973 IF ( pw%in_use /= COMPLEXDATA1D ) THEN 00974 CALL stop_program(routineN,moduleN,__LINE__,"Data field has to be COMPLEXDATA1D") 00975 ENDIF 00976 00977 IF ( pw%pw_grid%para%mode /= PW_MODE_DISTRIBUTED ) THEN 00978 CALL stop_program(routineN,moduleN,__LINE__,"This grid type is not distributed") 00979 ENDIF 00980 00981 ! after the gather we are in g-space 00982 pw%in_space = RECIPROCALSPACE 00983 00984 mapl => pw%pw_grid%mapl%pos 00985 mapm => pw%pw_grid%mapm%pos 00986 mapn => pw%pw_grid%mapn%pos 00987 00988 ngpts = SIZE ( pw%pw_grid%gsq ) 00989 00990 ghat => pw%pw_grid%g_hat 00991 yzq => pw%pw_grid%para%yzq 00992 00993 IF ( PRESENT ( scale ) ) THEN 00994 !$omp parallel do default(none), & 00995 !$omp private(l,m,n,mn), & 00996 !$omp shared(ngpts,mapl,mapm,mapn,ghat,yzq,pw,c,scale) 00997 DO gpt = 1, ngpts 00998 00999 l = mapl ( ghat ( 1, gpt ) ) + 1 01000 m = mapm ( ghat ( 2, gpt ) ) + 1 01001 n = mapn ( ghat ( 3, gpt ) ) + 1 01002 mn = yzq ( m, n ) 01003 pw%cc ( gpt ) = scale * c ( l, mn ) 01004 01005 END DO 01006 !$omp end parallel do 01007 ELSE 01008 !$omp parallel do default(none), & 01009 !$omp private(l,m,n,mn), & 01010 !$omp shared(ngpts,mapl,mapm,mapn,ghat,yzq,pw,c) 01011 DO gpt = 1, ngpts 01012 01013 l = mapl ( ghat ( 1, gpt ) ) + 1 01014 m = mapm ( ghat ( 2, gpt ) ) + 1 01015 n = mapn ( ghat ( 3, gpt ) ) + 1 01016 mn = yzq ( m, n ) 01017 pw%cc ( gpt ) = c ( l, mn ) 01018 01019 END DO 01020 !$omp end parallel do 01021 END IF 01022 01023 cpy = REAL ( ngpts,KIND=dp) * 1.e-6_dp 01024 CALL timestop(handle) 01025 01026 END SUBROUTINE pw_gather_p 01027 01028 ! ***************************************************************************** 01034 SUBROUTINE pw_scatter_s ( pw, c, scale, error) 01035 01036 TYPE(pw_type), INTENT(IN) :: pw 01037 COMPLEX(KIND=dp), DIMENSION(:, :, :), 01038 INTENT(INOUT) :: c 01039 REAL(KIND=dp), INTENT(IN), OPTIONAL :: scale 01040 TYPE(cp_error_type), INTENT(inout) :: error 01041 01042 CHARACTER(len=*), PARAMETER :: routineN = 'pw_scatter_s', 01043 routineP = moduleN//':'//routineN 01044 01045 INTEGER :: gpt, handle, l, m, n, ngpts 01046 INTEGER, DIMENSION(:), POINTER :: mapl, mapm, mapn 01047 INTEGER, DIMENSION(:, :), POINTER :: ghat 01048 LOGICAL :: failure 01049 REAL(KIND=dp) :: cpy 01050 01051 failure = .FALSE. 01052 CALL timeset(routineN,handle) 01053 CPPrecondition(pw%ref_count>0,cp_failure_level,routineP,error,failure) 01054 01055 IF ( pw%in_use /= COMPLEXDATA1D ) THEN 01056 CALL stop_program(routineN,moduleN,__LINE__,"Data field has to be COMPLEXDATA1D") 01057 ENDIF 01058 01059 IF ( pw%in_space /= RECIPROCALSPACE ) THEN 01060 CALL stop_program(routineN,moduleN,__LINE__,"Data has to be in RECIPROCALSPACE") 01061 ENDIF 01062 01063 mapl => pw%pw_grid%mapl%pos 01064 mapm => pw%pw_grid%mapm%pos 01065 mapn => pw%pw_grid%mapn%pos 01066 01067 ghat => pw%pw_grid%g_hat 01068 01069 ngpts = SIZE ( pw%pw_grid%gsq ) 01070 01071 ! should only zero the unused bits (but the zero is needed) 01072 IF ( .NOT. PRESENT ( scale ) ) c = 0.0_dp 01073 01074 IF ( PRESENT ( scale ) ) THEN 01075 01076 !$omp parallel do private(gpt,l,m,n) 01077 DO gpt = 1, ngpts 01078 01079 l = mapl ( ghat ( 1, gpt ) ) + 1 01080 m = mapm ( ghat ( 2, gpt ) ) + 1 01081 n = mapn ( ghat ( 3, gpt ) ) + 1 01082 c ( l, m, n ) = scale * pw%cc ( gpt ) 01083 01084 END DO 01085 01086 ELSE 01087 01088 !$omp parallel do private(gpt,l,m,n) 01089 DO gpt = 1, ngpts 01090 01091 l = mapl ( ghat ( 1, gpt ) ) + 1 01092 m = mapm ( ghat ( 2, gpt ) ) + 1 01093 n = mapn ( ghat ( 3, gpt ) ) + 1 01094 c ( l, m, n ) = pw%cc ( gpt ) 01095 01096 END DO 01097 01098 END IF 01099 01100 IF ( pw%pw_grid%grid_span == HALFSPACE ) THEN 01101 01102 mapl => pw%pw_grid%mapl%neg 01103 mapm => pw%pw_grid%mapm%neg 01104 mapn => pw%pw_grid%mapn%neg 01105 01106 IF ( PRESENT ( scale ) ) THEN 01107 01108 !$omp parallel do private(gpt,l,m,n) 01109 DO gpt = 1, ngpts 01110 01111 l = mapl ( ghat ( 1, gpt ) ) + 1 01112 m = mapm ( ghat ( 2, gpt ) ) + 1 01113 n = mapn ( ghat ( 3, gpt ) ) + 1 01114 c ( l, m, n ) = scale * CONJG ( pw%cc ( gpt ) ) 01115 01116 END DO 01117 01118 ELSE 01119 01120 !$omp parallel do private(gpt,l,m,n) 01121 DO gpt = 1, ngpts 01122 01123 l = mapl ( ghat ( 1, gpt ) ) + 1 01124 m = mapm ( ghat ( 2, gpt ) ) + 1 01125 n = mapn ( ghat ( 3, gpt ) ) + 1 01126 c ( l, m, n ) = CONJG ( pw%cc ( gpt ) ) 01127 01128 END DO 01129 01130 END IF 01131 01132 END IF 01133 01134 cpy = REAL ( ngpts,KIND=dp) * 1.e-6_dp 01135 CALL timestop(handle) 01136 01137 END SUBROUTINE pw_scatter_s 01138 01139 ! ***************************************************************************** 01140 SUBROUTINE pw_scatter_p ( pw, c, scale, error) 01141 01142 TYPE(pw_type), INTENT(IN), TARGET :: pw 01143 COMPLEX(KIND=dp), DIMENSION(:, :), 01144 INTENT(INOUT) :: c 01145 REAL(KIND=dp), INTENT(IN), OPTIONAL :: scale 01146 TYPE(cp_error_type), INTENT(inout) :: error 01147 01148 CHARACTER(len=*), PARAMETER :: routineN = 'pw_scatter_p', 01149 routineP = moduleN//':'//routineN 01150 01151 INTEGER :: gpt, handle, l, lb, m, mn, 01152 my_id, n, ngpts, num_threads, 01153 ub 01154 INTEGER, DIMENSION(:), POINTER :: mapl, mapm, mapn 01155 INTEGER, DIMENSION(:, :), POINTER :: ghat, yzq 01156 LOGICAL :: failure 01157 REAL(KIND=dp) :: cpy 01158 01159 !$ INTEGER :: omp_get_max_threads, omp_get_thread_num 01160 01161 failure = .FALSE. 01162 CALL timeset(routineN,handle) 01163 CPPrecondition(pw%ref_count>0,cp_failure_level,routineP,error,failure) 01164 01165 IF ( pw%in_use /= COMPLEXDATA1D ) THEN 01166 CALL stop_program(routineN,moduleN,__LINE__,"Data field has to be COMPLEXDATA1D") 01167 ENDIF 01168 01169 IF ( pw%in_space /= RECIPROCALSPACE ) THEN 01170 CALL stop_program(routineN,moduleN,__LINE__,"Data has to be in RECIPROCALSPACE") 01171 ENDIF 01172 01173 IF ( pw%pw_grid%para%mode /= PW_MODE_DISTRIBUTED ) THEN 01174 CALL stop_program(routineN,moduleN,__LINE__,"This grid type is not distributed") 01175 ENDIF 01176 01177 mapl => pw%pw_grid%mapl%pos 01178 mapm => pw%pw_grid%mapm%pos 01179 mapn => pw%pw_grid%mapn%pos 01180 01181 ghat => pw%pw_grid%g_hat 01182 yzq => pw%pw_grid%para%yzq 01183 01184 ngpts = SIZE ( pw%pw_grid%gsq ) 01185 01186 IF ( .NOT. PRESENT ( scale ) ) THEN 01187 num_threads = 1 01188 my_id = 0 01189 !$omp parallel default(none), & 01190 !$omp private(my_id, num_threads, lb, ub), & 01191 !$omp shared(c) 01192 !$ num_threads = MIN(omp_get_max_threads(), SIZE(c,2)) 01193 !$ my_id = omp_get_thread_num() 01194 IF (my_id < num_threads) THEN 01195 lb = (SIZE(c,2)*my_id)/num_threads + 1 01196 ub = (SIZE(c,2)*(my_id+1))/num_threads 01197 c(:, lb:ub) = 0.0_dp 01198 END IF 01199 !$omp end parallel 01200 END IF 01201 01202 IF ( PRESENT ( scale ) ) THEN 01203 !$omp parallel do default(none), & 01204 !$omp private(l,m,n,mn), & 01205 !$omp shared(ngpts,mapl,mapm,mapn,ghat,yzq,pw,c,scale) 01206 DO gpt = 1, ngpts 01207 01208 l = mapl ( ghat ( 1, gpt ) ) + 1 01209 m = mapm ( ghat ( 2, gpt ) ) + 1 01210 n = mapn ( ghat ( 3, gpt ) ) + 1 01211 mn = yzq ( m, n ) 01212 c ( l, mn ) = scale * pw%cc ( gpt ) 01213 01214 END DO 01215 !$omp end parallel do 01216 ELSE 01217 !$omp parallel do default(none), & 01218 !$omp private(l,m,n,mn), & 01219 !$omp shared(ngpts,mapl,mapm,mapn,ghat,yzq,pw,c) 01220 DO gpt = 1, ngpts 01221 01222 l = mapl ( ghat ( 1, gpt ) ) + 1 01223 m = mapm ( ghat ( 2, gpt ) ) + 1 01224 n = mapn ( ghat ( 3, gpt ) ) + 1 01225 mn = yzq ( m, n ) 01226 c ( l, mn ) = pw%cc ( gpt ) 01227 01228 END DO 01229 !$omp end parallel do 01230 END IF 01231 01232 IF ( pw%pw_grid%grid_span == HALFSPACE ) THEN 01233 01234 mapm => pw%pw_grid%mapm%neg 01235 mapn => pw%pw_grid%mapn%neg 01236 mapl => pw%pw_grid%mapl%neg 01237 01238 IF ( PRESENT ( scale ) ) THEN 01239 !$omp parallel do default(none), & 01240 !$omp private(l,m,n,mn), & 01241 !$omp shared(ngpts,mapl,mapm,mapn,ghat,yzq,pw,c,scale) 01242 DO gpt = 1, ngpts 01243 01244 l = mapl ( ghat ( 1, gpt ) ) + 1 01245 m = mapm ( ghat ( 2, gpt ) ) + 1 01246 n = mapn ( ghat ( 3, gpt ) ) + 1 01247 mn = yzq ( m, n ) 01248 c ( l, mn ) = scale * CONJG ( pw%cc ( gpt ) ) 01249 01250 END DO 01251 !$omp end parallel do 01252 ELSE 01253 !$omp parallel do default(none), & 01254 !$omp private(l,m,n,mn), & 01255 !$omp shared(ngpts,mapl,mapm,mapn,ghat,yzq,pw,c) 01256 DO gpt = 1, ngpts 01257 01258 l = mapl ( ghat ( 1, gpt ) ) + 1 01259 m = mapm ( ghat ( 2, gpt ) ) + 1 01260 n = mapn ( ghat ( 3, gpt ) ) + 1 01261 mn = yzq ( m, n ) 01262 c ( l, mn ) = CONJG ( pw%cc ( gpt ) ) 01263 01264 END DO 01265 !$omp end parallel do 01266 END IF 01267 01268 END IF 01269 01270 cpy = REAL ( ngpts,KIND=dp) * 1.e-6_dp 01271 CALL timestop(handle) 01272 01273 01274 END SUBROUTINE pw_scatter_p 01275 01276 ! ***************************************************************************** 01287 SUBROUTINE fft_wrap_pw1 ( pw1, debug, error) 01288 01289 TYPE(pw_type), INTENT(INOUT), TARGET :: pw1 01290 LOGICAL, INTENT(IN), OPTIONAL :: debug 01291 TYPE(cp_error_type), INTENT(inout) :: error 01292 01293 CHARACTER(len=*), PARAMETER :: routineN = 'fft_wrap_pw1', 01294 routineP = moduleN//':'//routineN 01295 01296 INTEGER :: dir, handle, out_space 01297 LOGICAL :: failure, test 01298 REAL(KIND=dp) :: norm 01299 01300 failure = .FALSE. 01301 CALL timeset(routineN,handle) 01302 CPPrecondition(pw1%ref_count>0,cp_failure_level,routineP,error,failure) 01303 01304 IF ( PRESENT ( debug ) ) THEN 01305 test = debug 01306 ELSE 01307 test = .FALSE. 01308 END IF 01309 01310 IF ( pw1%in_use == COMPLEXDATA3D .AND. & 01311 pw1%pw_grid%para%mode == PW_MODE_LOCAL ) THEN 01312 !..what dirction will the transform be? 01313 IF ( pw1%in_space == REALSPACE ) THEN 01314 dir = FWFFT 01315 norm = pw1%pw_grid%dvol 01316 out_space = RECIPROCALSPACE 01317 ELSEIF ( pw1%in_space == RECIPROCALSPACE ) THEN 01318 dir = BWFFT 01319 norm = 1.0_dp 01320 out_space = REALSPACE 01321 ELSE 01322 CALL stop_program(routineN,moduleN,__LINE__, "PW structure is missing a "//& 01323 "proper tag to identidy its space" ) 01324 END IF 01325 01326 CALL fft3d ( dir, pw1%pw_grid%npts, pw1%cc3d, & 01327 scale = norm, debug=test ) 01328 01329 !..tag new data with correct space 01330 pw1%in_space = out_space 01331 ELSE 01332 CALL stop_program(routineN,moduleN,__LINE__, "In place FFT only possible for "//& 01333 "replicated data with COMPLEXDATA3D structure" ) 01334 END IF 01335 CALL timestop(handle) 01336 01337 END SUBROUTINE fft_wrap_pw1 01338 01339 ! ***************************************************************************** 01340 SUBROUTINE fft_wrap_pw1pw2 ( pw1, pw2, debug, error) 01341 01342 TYPE(pw_type), INTENT(IN), TARGET :: pw1 01343 TYPE(pw_type), INTENT(INOUT), TARGET :: pw2 01344 LOGICAL, INTENT(IN), OPTIONAL :: debug 01345 TYPE(cp_error_type), INTENT(inout) :: error 01346 01347 CHARACTER(len=*), PARAMETER :: routineN = 'fft_wrap_pw1pw2', 01348 routineP = moduleN//':'//routineN 01349 01350 CHARACTER(LEN=9) :: mode 01351 COMPLEX(KIND=dp), DIMENSION(:, :), 01352 POINTER :: grays 01353 COMPLEX(KIND=dp), DIMENSION(:, :, :), 01354 POINTER :: c_in, c_out 01355 INTEGER :: dir, handle, handle2, l1, l2, l3, lb, my_id, my_pos, 01356 nloc( 3 ), nrays, nsize, num_threads, out_space, stat, ub 01357 INTEGER, DIMENSION(:), POINTER :: n 01358 LOGICAL :: failure, test 01359 REAL(KIND=dp) :: norm 01360 01361 !$ INTEGER :: omp_get_max_threads, omp_get_thread_num 01362 01363 failure = .FALSE. 01364 CALL timeset(routineN,handle2) 01365 CALL timeset(routineN//"_"//TRIM(ADJUSTL(cp_to_string( & 01366 CEILING (pw1%pw_grid%cutoff/10)*10))),handle) 01367 01368 CPPrecondition(pw1%ref_count>0,cp_failure_level,routineP,error,failure) 01369 CPPrecondition(pw2%ref_count>0,cp_failure_level,routineP,error,failure) 01370 NULLIFY ( c_in ) 01371 NULLIFY ( c_out ) 01372 01373 IF ( PRESENT ( debug ) ) THEN 01374 test = debug 01375 ELSE 01376 test = .FALSE. 01377 END IF 01378 01379 !..check if grids are compatible 01380 IF ( pw1%pw_grid %id_nr /= pw2%pw_grid %id_nr ) THEN 01381 IF ( pw1%pw_grid%dvol /= pw2%pw_grid%dvol ) THEN 01382 CALL stop_program(routineN,moduleN,__LINE__,& 01383 "PW grids not compatible") 01384 END IF 01385 IF ( pw1%pw_grid%para %group /= pw2%pw_grid%para%group) THEN 01386 CALL stop_program(routineN,moduleN,__LINE__,& 01387 "PW grids have not compatible MPI groups") 01388 END IF 01389 END IF 01390 01391 !..prepare input 01392 IF ( pw1%in_space == REALSPACE ) THEN 01393 dir = FWFFT 01394 norm = 1.0_dp / pw1%pw_grid%ngpts 01395 out_space = RECIPROCALSPACE 01396 ELSE IF ( pw1%in_space == RECIPROCALSPACE ) THEN 01397 dir = BWFFT 01398 norm = 1.0_dp 01399 out_space = REALSPACE 01400 ELSE 01401 CALL stop_program(routineN,moduleN,__LINE__,"Error in space tag") 01402 END IF 01403 01404 n => pw1%pw_grid%npts 01405 01406 mode = fftselect ( pw1%in_use, pw2%in_use, pw1%in_space, error=error) 01407 01408 IF ( pw1%pw_grid%para%mode == PW_MODE_LOCAL ) THEN 01409 01410 ! 01411 !..replicated data, use local FFT 01412 ! 01413 01414 IF ( test ) THEN 01415 WRITE ( *,'(A)') " FFT Protocol " 01416 IF ( dir == FWFFT ) WRITE ( *,'(A,T76,A)') " Transform direction ","FWFFT" 01417 IF ( dir == BWFFT ) WRITE ( *,'(A,T76,A)') " Transform direction ","BWFFT" 01418 IF ( pw1%in_space == REALSPACE ) & 01419 WRITE ( *,'(A,T72,A)') " in space ","REALSPACE" 01420 IF ( pw1%in_space == RECIPROCALSPACE ) & 01421 WRITE ( *,'(A,T66,A)') " in space ","RECIPROCALSPACE" 01422 IF ( out_space == REALSPACE ) & 01423 WRITE ( *,'(A,T72,A)') " out space ","REALSPACE" 01424 IF ( out_space == RECIPROCALSPACE ) & 01425 WRITE ( *,'(A,T66,A)') " out space ","RECIPROCALSPACE" 01426 WRITE ( *,'(A,T66,E15.6)') " scale factor",norm 01427 END IF 01428 01429 SELECT CASE ( mode ) 01430 CASE DEFAULT 01431 CALL stop_program(routineN,moduleN,__LINE__,& 01432 "Illegal combination of in_use and in_space") 01433 CASE ( "FW_C3DC3D" ) 01434 c_in => pw1%cc3d 01435 c_out => pw2%cc3d 01436 CALL fft3d ( dir, n, c_in, c_out, scale = norm, debug=test ) 01437 CASE ( "FW_R3DC3D" ) 01438 pw2%cc3d = CMPLX ( pw1%cr3d, 0.0_dp, KIND=dp) 01439 c_out => pw2%cc3d 01440 CALL fft3d ( dir, n, c_out, scale = norm, debug=test ) 01441 CASE ( "FW_C3DC1D" ) 01442 c_in => pw1%cc3d 01443 ALLOCATE ( c_out( n(1), n(2), n(3) ), STAT = stat) 01444 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01445 ! transform 01446 CALL fft3d ( dir, n, c_in, c_out, scale = norm, debug=test ) 01447 ! gather results 01448 IF ( test ) WRITE ( *,'(A)') " PW_GATHER : 3d -> 1d " 01449 CALL pw_gather ( pw2, c_out, error=error) 01450 DEALLOCATE ( c_out, STAT = stat ) 01451 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01452 CASE ( "FW_R3DC1D" ) 01453 #if defined ( __CUDAPW ) 01454 CALL cuda_pw_fft_wrap_r3dc1d(pw1, pw2, dir, n, norm) 01455 #else 01456 ALLOCATE ( c_out( n(1), n(2), n(3) ), STAT = stat ) 01457 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01458 nsize=SIZE(pw1%cr3d,1)*SIZE(pw1%cr3d,2)*SIZE(pw1%cr3d,3) 01459 l1 = LBOUND(pw1%cr3d,1) 01460 l2 = LBOUND(pw1%cr3d,2) 01461 l3 = LBOUND(pw1%cr3d,3) 01462 CALL copy_rc(nsize,pw1%cr3d(l1,l2,l3),c_out(1,1,1)) 01463 CALL fft3d ( dir, n, c_out, scale = 1._dp, debug=test ) 01464 CALL pw_gather ( pw2, c_out, scale = norm, error=error) 01465 DEALLOCATE ( c_out, STAT = stat ) 01466 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01467 #endif 01468 CASE ( "BW_C3DC3D" ) 01469 c_in => pw1%cc3d 01470 c_out => pw2%cc3d 01471 CALL fft3d ( dir, n, c_in, c_out, scale = norm, debug=test ) 01472 CASE ( "BW_C3DR3D" ) 01473 c_in => pw1%cc3d 01474 ALLOCATE ( c_out( n(1), n(2), n(3) ), STAT = stat ) 01475 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01476 CALL fft3d ( dir, n, c_in, c_out, scale = norm, debug=test ) 01477 ! use real part only 01478 IF ( test ) WRITE ( *,'(A)') " REAL part " 01479 pw2%cr3d = REAL ( c_out,KIND=dp) 01480 DEALLOCATE ( c_out, STAT = stat ) 01481 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01482 CASE ( "BW_C1DC3D" ) 01483 c_out => pw2%cc3d 01484 IF ( test ) WRITE ( *,'(A)') " PW_SCATTER : 3d -> 1d " 01485 CALL pw_scatter ( pw1, c_out, error=error) 01486 CALL fft3d ( dir, n, c_out, scale = norm, debug=test ) 01487 CASE ( "BW_C1DR3D" ) 01488 #if defined ( __CUDAPW ) 01489 CALL cuda_pw_fft_wrap_c1dr3d(pw1, pw2, dir, n, norm) 01490 #else 01491 ALLOCATE ( c_out( n(1), n(2), n(3) ), STAT = stat) 01492 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01493 IF ( test ) WRITE ( *,'(A)') " PW_SCATTER : 3d -> 1d " 01494 CALL pw_scatter ( pw1, c_out, error=error) 01495 ! transform 01496 CALL fft3d ( dir, n, c_out, scale = norm, debug=test ) 01497 ! use real part only 01498 IF ( test ) WRITE ( *,'(A)') " REAL part " 01499 ! pw2%cr3d = REAL ( c_out,KIND=dp) 01500 nsize=SIZE(pw2%cr3d,1)*SIZE(pw2%cr3d,2)*SIZE(pw2%cr3d,3) 01501 l1 = LBOUND(pw2%cr3d,1) 01502 l2 = LBOUND(pw2%cr3d,2) 01503 l3 = LBOUND(pw2%cr3d,3) 01504 CALL copy_cr(nsize,c_out(1,1,1),pw2%cr3d(l1,l2,l3)) 01505 DEALLOCATE ( c_out, STAT = stat ) 01506 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01507 #endif 01508 END SELECT 01509 01510 IF ( test ) WRITE ( *,'(A)') " End of FFT Protocol " 01511 01512 ELSE 01513 01514 ! 01515 !..parallel FFT 01516 ! 01517 01518 IF ( test .AND. pw1%pw_grid%para%group_head ) THEN 01519 WRITE ( *,'(A)') " FFT Protocol " 01520 IF ( dir == FWFFT ) WRITE ( *,'(A,T76,A)') " Transform direction ","FWFFT" 01521 IF ( pw1%in_space == REALSPACE ) & 01522 WRITE ( *,'(A,T72,A)') " in space ","REALSPACE" 01523 IF ( pw1%in_space == RECIPROCALSPACE ) & 01524 WRITE ( *,'(A,T66,A)') " in space ","RECIPROCALSPACE" 01525 IF ( out_space == REALSPACE ) & 01526 WRITE ( *,'(A,T72,A)') " out space ","REALSPACE" 01527 IF ( out_space == RECIPROCALSPACE ) & 01528 WRITE ( *,'(A,T66,A)') " out space ","RECIPROCALSPACE" 01529 WRITE ( *,'(A,T66,E15.6)') " scale factor",norm 01530 END IF 01531 01532 my_pos = pw1%pw_grid%para%my_pos 01533 nrays = pw1%pw_grid%para%nyzray ( my_pos ) 01534 grays => pw1%pw_grid%grays 01535 CPPostcondition(SIZE(grays,1)==n(1),cp_failure_level,routineP,error,failure) 01536 CPPostcondition(SIZE(grays,2)==nrays,cp_failure_level,routineP,error,failure) 01537 01538 num_threads = 1 01539 my_id = 0 01540 !$omp parallel default(none), & 01541 !$omp private(my_id, num_threads, lb, ub), & 01542 !$omp shared(grays,nrays) 01543 !$ num_threads = MIN(omp_get_max_threads(), nrays) 01544 !$ my_id = omp_get_thread_num() 01545 IF (my_id < num_threads) THEN 01546 lb = (nrays*my_id)/num_threads + 1 01547 ub = (nrays*(my_id+1))/num_threads 01548 grays(:, lb:ub) = 0.0_dp 01549 END IF 01550 !$omp end parallel 01551 01552 SELECT CASE ( mode ) 01553 CASE DEFAULT 01554 CALL stop_program(routineN,moduleN,__LINE__,& 01555 "Illegal combination of in_use and in_space "//& 01556 "in parallel 3d FFT") 01557 CASE ( "FW_C3DC1D" ) 01558 c_in => pw1%cc3d 01559 CASE ( "FW_R3DC1D" ) 01560 nloc = pw1%pw_grid%npts_local 01561 ALLOCATE ( c_in( nloc(1), nloc(2), nloc(3) ), STAT = stat) 01562 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01563 c_in = CMPLX ( pw1%cr3d, 0.0_dp,KIND=dp) 01564 CASE ( "BW_C1DC3D" ) 01565 IF ( test .AND. pw1%pw_grid%para%group_head ) & 01566 WRITE ( *,'(A)') " PW_SCATTER : 2d -> 1d " 01567 CALL pw_scatter ( pw1, grays, error=error) 01568 c_in => pw2%cc3d 01569 CASE ( "BW_C1DR3D" ) 01570 IF ( test .AND. pw1%pw_grid%para%group_head ) & 01571 WRITE ( *,'(A)') " PW_SCATTER : 2d -> 1d " 01572 CALL pw_scatter ( pw1, grays, error=error) 01573 nloc = pw2%pw_grid%npts_local 01574 ALLOCATE ( c_in( nloc(1), nloc(2), nloc(3) ), STAT = stat ) 01575 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01576 END SELECT 01577 01578 !..transform 01579 IF ( pw1%pw_grid%para%ray_distribution ) THEN 01580 CALL fft3d ( dir, n, c_in, grays, pw1%pw_grid%para%group, & 01581 pw1%pw_grid%para%rs_group, & 01582 pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, & 01583 pw1%pw_grid%para%bo, scale = norm, debug=test ) 01584 ELSE 01585 CALL fft3d ( dir, n, c_in, grays, pw1%pw_grid%para%rs_group, & 01586 pw1%pw_grid%para%bo, scale = norm, debug=test ) 01587 END IF 01588 01589 !..prepare output 01590 SELECT CASE ( mode ) 01591 CASE DEFAULT 01592 CALL stop_program(routineN,moduleN,__LINE__, & 01593 "Illegal combination of in_use and in_space "//& 01594 "in parallel 3d FFT" ) 01595 CASE ( "FW_C3DC1D" ) 01596 IF ( test .AND. pw1%pw_grid%para%group_head ) & 01597 WRITE ( *,'(A)') " PW_GATHER : 2d -> 1d " 01598 CALL pw_gather ( pw2, grays, error=error) 01599 CASE ( "FW_R3DC1D" ) 01600 IF ( test .AND. pw1%pw_grid%para%group_head ) & 01601 WRITE ( *,'(A)') " PW_GATHER : 2d -> 1d " 01602 CALL pw_gather ( pw2, grays, error=error) 01603 DEALLOCATE ( c_in, STAT = stat ) 01604 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01605 CASE ( "BW_C1DC3D" ) 01606 ! nothing to do 01607 CASE ( "BW_C1DR3D" ) 01608 IF ( test .AND. pw1%pw_grid%para%group_head ) & 01609 WRITE ( *,'(A)') " Real part " 01610 pw2%cr3d = REAL ( c_in,KIND=dp) 01611 DEALLOCATE ( c_in, STAT = stat ) 01612 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01613 END SELECT 01614 01615 END IF 01616 01617 !..update the space tag for pw2 01618 pw2%in_space = out_space 01619 01620 IF ( test .AND. pw1%pw_grid%para%group_head ) THEN 01621 WRITE ( *,'(A)') " End of FFT Protocol " 01622 END IF 01623 CALL timestop(handle) 01624 CALL timestop(handle2) 01625 01626 01627 END SUBROUTINE fft_wrap_pw1pw2 01628 01629 ! ***************************************************************************** 01630 FUNCTION fftselect ( use1, use2, space1, error) RESULT ( mode ) 01631 INTEGER, INTENT(IN) :: use1, use2, space1 01632 TYPE(cp_error_type), INTENT(inout) :: error 01633 CHARACTER(LEN=9) :: mode 01634 01635 CHARACTER(len=*), PARAMETER :: routineN = 'fftselect', 01636 routineP = moduleN//':'//routineN 01637 01638 LOGICAL :: failure 01639 01640 failure = .FALSE. 01641 IF ( space1 == REALSPACE ) THEN 01642 mode ( 1 : 3 ) = "FW_" 01643 ELSE IF ( space1 == RECIPROCALSPACE ) THEN 01644 mode ( 1 : 3 ) = "BW_" 01645 ELSE 01646 CALL stop_program(routineN,moduleN,__LINE__,"Error in space tag") 01647 END IF 01648 01649 SELECT CASE ( use1 ) 01650 CASE ( COMPLEXDATA3D ) 01651 mode ( 4 : 6 ) = "C3D" 01652 CASE ( REALDATA3D ) 01653 mode ( 4 : 6 ) = "R3D" 01654 CASE ( COMPLEXDATA1D ) 01655 mode ( 4 : 6 ) = "C1D" 01656 CASE ( REALDATA1D ) 01657 mode ( 4 : 6 ) = "R1D" 01658 CASE DEFAULT 01659 CALL stop_program(routineN,moduleN,__LINE__,"Error in use1 tag") 01660 END SELECT 01661 01662 SELECT CASE ( use2 ) 01663 CASE ( COMPLEXDATA3D ) 01664 mode ( 7 : 9 ) = "C3D" 01665 CASE ( REALDATA3D ) 01666 mode ( 7 : 9 ) = "R3D" 01667 CASE ( COMPLEXDATA1D ) 01668 mode ( 7 : 9 ) = "C1D" 01669 CASE ( REALDATA1D ) 01670 mode ( 7 : 9 ) = "R1D" 01671 CASE DEFAULT 01672 CALL stop_program(routineN,moduleN,__LINE__,"Error in use1 tag") 01673 END SELECT 01674 01675 END FUNCTION fftselect 01676 01677 ! ***************************************************************************** 01689 SUBROUTINE pw_write(pw, unit_nr, error) 01690 TYPE(pw_type), INTENT(in) :: pw 01691 INTEGER, INTENT(in) :: unit_nr 01692 TYPE(cp_error_type), INTENT(inout) :: error 01693 01694 CHARACTER(len=*), PARAMETER :: routineN = 'pw_write', 01695 routineP = moduleN//':'//routineN 01696 01697 INTEGER :: iostat 01698 LOGICAL :: failure 01699 01700 failure=.FALSE. 01701 01702 WRITE (unit=unit_nr, fmt="('<pw>:{ id_nr=',i8,',')",iostat=iostat)& 01703 pw%id_nr 01704 01705 SELECT CASE(pw%in_use) 01706 CASE (REALDATA1D) 01707 WRITE (unit=unit_nr, fmt="(a)",iostat=iostat) " in_use=REALDATA1D" 01708 IF (ASSOCIATED(pw%cr)) THEN 01709 WRITE (unit=unit_nr, fmt="(' cr=<real(',i8,':',i8,')at 0x',z16.16,'>')")& 01710 LBOUND(pw%cr,1),UBOUND(pw%cr,1),m_loc_r(pw%cr(LBOUND(pw%cr))) 01711 ELSE 01712 WRITE (unit=unit_nr, fmt="(' cr=*null*')") 01713 END IF 01714 CASE (REALDATA3D) 01715 WRITE (unit=unit_nr, fmt="(a)",iostat=iostat) " in_use=REALDATA3D" 01716 IF (ASSOCIATED(pw%cr3d)) THEN 01717 WRITE (unit=unit_nr, fmt="(' cr3d=<real(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')at 0x',z16.16,'>')")& 01718 LBOUND(pw%cr3d,1),UBOUND(pw%cr3d,1),LBOUND(pw%cr3d,2),UBOUND(pw%cr3d,2),& 01719 LBOUND(pw%cr3d,3),UBOUND(pw%cr3d,3),& 01720 m_loc_r(pw%cr3d) 01721 ELSE 01722 WRITE (unit=unit_nr, fmt="(' cr3d=*null*')") 01723 END IF 01724 CASE (COMPLEXDATA1D) 01725 WRITE (unit=unit_nr, fmt="(a)",iostat=iostat) " in_use=COMPLEXDATA1D" 01726 IF (ASSOCIATED(pw%cc)) THEN 01727 WRITE (unit=unit_nr, fmt="(' cc=<real(',i8,':',i8,') at 0x',z16.16,'>')")& 01728 LBOUND(pw%cc,1),UBOUND(pw%cc,1),m_loc_c(pw%cc) 01729 ELSE 01730 WRITE (unit=unit_nr, fmt="(' cc=*null*')") 01731 END IF 01732 CASE (COMPLEXDATA3D) 01733 WRITE (unit=unit_nr, fmt="(a)",iostat=iostat) " in_use=COMPLEXDATA3D" 01734 IF (ASSOCIATED(pw%cc3d)) THEN 01735 WRITE (unit=unit_nr, fmt="(' cc3d=<real(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,') at 0x',z16.16,'>')")& 01736 LBOUND(pw%cc3d,1),UBOUND(pw%cc3d,1),LBOUND(pw%cc3d,2),UBOUND(pw%cc3d,2),& 01737 LBOUND(pw%cc3d,3),UBOUND(pw%cc3d,3),& 01738 m_loc_c(pw%cc3d) 01739 ELSE 01740 WRITE (unit=unit_nr, fmt="(' cr3d=*null*')") 01741 END IF 01742 CASE default 01743 WRITE (unit=unit_nr, fmt="(' in_use=',i8,',')",iostat=iostat)& 01744 pw%in_use 01745 END SELECT 01746 01747 SELECT CASE(pw%in_space) 01748 CASE (NOSPACE) 01749 WRITE (unit=unit_nr, fmt="(a)",iostat=iostat) " in_space=NOSPACE" 01750 CASE (REALSPACE) 01751 WRITE (unit=unit_nr, fmt="(a)",iostat=iostat) " in_space=REALSPACE" 01752 CASE (RECIPROCALSPACE) 01753 WRITE (unit=unit_nr, fmt="(a)",iostat=iostat) " in_space=RECIPROCALSPACE" 01754 CASE default 01755 WRITE (unit=unit_nr, fmt="(' in_space=',i8,',')",iostat=iostat)& 01756 pw%in_space 01757 END SELECT 01758 01759 WRITE (unit=unit_nr, fmt="(' pw_grid%id_nr=',i8,/,' }')",iostat=iostat)& 01760 pw%pw_grid%id_nr 01761 01762 END SUBROUTINE pw_write 01763 01764 ! ***************************************************************************** 01765 FUNCTION pw_compatible ( grida, gridb, error) RESULT ( compat ) 01766 TYPE(pw_grid_type), INTENT(IN) :: grida, gridb 01767 TYPE(cp_error_type), INTENT(inout) :: error 01768 LOGICAL :: compat 01769 01770 CHARACTER(len=*), PARAMETER :: routineN = 'pw_compatible', 01771 routineP = moduleN//':'//routineN 01772 01773 LOGICAL :: failure 01774 01775 failure = .FALSE. 01776 compat = .FALSE. 01777 IF ( grida%id_nr == gridb%id_nr ) THEN 01778 compat = .TRUE. 01779 ELSE IF ( grida%reference == gridb%id_nr ) THEN 01780 compat = .TRUE. 01781 ELSE IF ( gridb%reference == grida%id_nr ) THEN 01782 compat = .TRUE. 01783 END IF 01784 01785 END FUNCTION pw_compatible 01786 01787 ! ***************************************************************************** 01796 SUBROUTINE pw_compare_debug(pw1,pw2,maxdiff,error) 01797 TYPE(pw_type), POINTER :: pw1, pw2 01798 REAL(KIND=dp), INTENT(out), OPTIONAL :: maxdiff 01799 TYPE(cp_error_type), INTENT(inout) :: error 01800 01801 CHARACTER(len=*), PARAMETER :: routineN = 'pw_compare_debug', 01802 routineP = moduleN//':'//routineN 01803 01804 INTEGER :: i, j, k, unit_nr 01805 LOGICAL :: failure 01806 REAL(KIND=dp) :: diff, mdiff 01807 TYPE(cp_logger_type), POINTER :: logger 01808 01809 failure=.FALSE. 01810 01811 CPPrecondition(ASSOCIATED(pw1),cp_failure_level,routineP,error,failure) 01812 CPPrecondition(ASSOCIATED(pw2),cp_failure_level,routineP,error,failure) 01813 CPPrecondition(pw1%ref_count>0,cp_failure_level,routineP,error,failure) 01814 CPPrecondition(pw2%ref_count>0,cp_failure_level,routineP,error,failure) 01815 CPPrecondition(pw1%in_space==pw2%in_space,cp_warning_level,routineP,error,failure) 01816 CPPrecondition(pw1%in_use==pw2%in_use,cp_warning_level,routineP,error,failure) 01817 CALL cp_assert(ALL(pw1%pw_grid%bounds_local==pw2%pw_grid%bounds_local),& 01818 cp_failure_level,cp_assertion_failed,routineP,& 01819 "wrong pw distribution",error,failure) 01820 IF (.NOT. failure) THEN 01821 logger => cp_error_get_logger(error) 01822 unit_nr=-1 01823 mdiff=0.0_dp 01824 SELECT CASE(pw1%in_use) 01825 CASE(REALDATA3D) 01826 DO k=pw1%pw_grid%bounds_local(1,3),pw1%pw_grid%bounds_local(2,3) 01827 DO j=pw1%pw_grid%bounds_local(1,2),pw1%pw_grid%bounds_local(2,2) 01828 DO i=pw1%pw_grid%bounds_local(1,1),pw1%pw_grid%bounds_local(2,1) 01829 diff=ABS(pw1%cr3d(i,j,k)-pw2%cr3d(i,j,k)) 01830 IF (mdiff<diff) THEN 01831 WRITE(unit=cp_logger_get_default_unit_nr(logger,local=.TRUE.),& 01832 fmt="(' diff=',e10.4,'at(',i4,i4,i4,'),',e10.4,'vs',e10.4)")& 01833 diff,i,j,k,pw1%cr3d(i,j,k),pw2%cr3d(i,j,k) 01834 mdiff=diff 01835 END IF 01836 END DO 01837 END DO 01838 END DO 01839 CASE(COMPLEXDATA1D) 01840 DO i=1,pw1%pw_grid%ngpts_local 01841 diff=ABS(pw1%cc(i)-pw2%cc(i)) 01842 IF (mdiff<diff) THEN 01843 WRITE(unit=cp_logger_get_default_unit_nr(logger,local=.TRUE.),& 01844 fmt="(' diff=',e10.4,'at(',i4,i4,i4,'),',e9.3,e9.3,'vs',e9.3,e9.3)")& 01845 diff,i,pw1%cc(i),pw2%cc(i) 01846 mdiff=diff 01847 END IF 01848 END DO 01849 CASE default 01850 CPPrecondition(.FALSE.,cp_warning_level,routineP,error,failure) 01851 END SELECT 01852 WRITE(unit=cp_logger_get_default_unit_nr(logger,local=.TRUE.),& 01853 fmt="(' maxdiff=',e10.4)") mdiff 01854 IF (PRESENT(maxdiff)) maxdiff=mdiff 01855 ELSE 01856 WRITE(unit=cp_logger_get_default_unit_nr(logger,local=.TRUE.),& 01857 fmt="(' incompatible pws')") 01858 IF (PRESENT(maxdiff)) maxdiff=HUGE(0.0_dp) 01859 END IF 01860 END SUBROUTINE pw_compare_debug 01861 01862 ! ***************************************************************************** 01869 FUNCTION pw_integral_ab ( pw1, pw2, error) RESULT ( integral_value ) 01870 01871 TYPE(pw_type), INTENT(IN) :: pw1, pw2 01872 TYPE(cp_error_type), INTENT(inout) :: error 01873 REAL(KIND=dp) :: integral_value 01874 01875 CHARACTER(len=*), PARAMETER :: routineN = 'pw_integral_ab', 01876 routineP = moduleN//':'//routineN 01877 01878 LOGICAL :: failure 01879 01880 failure = .FALSE. 01881 IF ( pw1%pw_grid %id_nr /= pw2%pw_grid %id_nr ) THEN 01882 CALL stop_program(routineN,moduleN,__LINE__,"Grids incompatible") 01883 END IF 01884 01885 ! since the return value is real, only do accurate sum on the real bit ? 01886 IF ( pw1%in_use == REALDATA3D .AND. pw2%in_use == REALDATA3D ) THEN 01887 integral_value = accurate_sum ( pw1%cr3d ( :, :, : ) & 01888 * pw2%cr3d ( :, :, : ) ) 01889 ELSE IF ( pw1%in_use == REALDATA3D & 01890 .AND. pw2%in_use == COMPLEXDATA3D ) THEN 01891 integral_value = accurate_sum ( pw1%cr3d ( :, :, : ) & 01892 * pw2%cc3d ( :, :, : ) ) 01893 ELSE IF ( pw1%in_use == COMPLEXDATA3D & 01894 .AND. pw2%in_use == REALDATA3D ) THEN 01895 integral_value = accurate_sum ( pw1%cc3d ( :, :, : ) & 01896 * pw2%cr3d ( :, :, : ) ) 01897 ELSE IF ( pw1%in_use == COMPLEXDATA3D & 01898 .AND. pw2%in_use == COMPLEXDATA3D ) THEN 01899 integral_value = accurate_sum ( CONJG ( pw1%cc3d ( :, :, : ) ) & 01900 * pw2%cc3d ( :, :, : ) ) 01901 01902 ELSE IF ( pw1%in_use == REALDATA1D & 01903 .AND. pw2%in_use == REALDATA1D ) THEN 01904 integral_value = accurate_sum ( pw1%cr (:) * pw2%cr (:) ) 01905 ELSE IF ( pw1%in_use == REALDATA1D & 01906 .AND. pw2%in_use == COMPLEXDATA1D ) THEN 01907 integral_value = accurate_sum ( pw1%cr (:) * pw2%cc (:) ) 01908 ELSE IF ( pw1%in_use == COMPLEXDATA1D & 01909 .AND. pw2%in_use == REALDATA1D ) THEN 01910 integral_value = accurate_sum ( pw1%cc (:) * pw2%cr (:) ) 01911 ELSE IF ( pw1%in_use == COMPLEXDATA1D & 01912 .AND. pw2%in_use == COMPLEXDATA1D ) THEN 01913 integral_value = accurate_sum ( CONJG ( pw1%cc (:) ) * pw2%cc (:) ) 01914 ELSE 01915 CALL stop_program(routineN,moduleN,__LINE__,"No possible DATA") 01916 END IF 01917 01918 IF ( pw1%in_use == REALDATA3D .OR. pw1%in_use == COMPLEXDATA3D ) THEN 01919 integral_value = integral_value * pw1%pw_grid%dvol 01920 ELSE 01921 integral_value = integral_value * pw1%pw_grid%vol 01922 ENDIF 01923 IF ( pw1%in_use == COMPLEXDATA1D ) THEN 01924 IF ( pw1%pw_grid%grid_span == HALFSPACE ) THEN 01925 integral_value = 2.0_dp * integral_value 01926 IF ( pw1%pw_grid%have_g0 ) integral_value = integral_value - & 01927 CONJG ( pw1%cc ( 1 ) ) * pw2%cc ( 1 ) 01928 END IF 01929 END IF 01930 01931 IF ( pw1%pw_grid%para%mode == PW_MODE_DISTRIBUTED ) & 01932 CALL mp_sum ( integral_value, pw1%pw_grid%para%group ) 01933 01934 END FUNCTION pw_integral_ab 01935 01936 ! ***************************************************************************** 01943 FUNCTION pw_integral_ab_fast ( pw1, pw2, error) RESULT ( integral_value ) 01944 01945 TYPE(pw_type), INTENT(IN) :: pw1, pw2 01946 TYPE(cp_error_type), INTENT(inout) :: error 01947 REAL(KIND=dp) :: integral_value 01948 01949 CHARACTER(len=*), PARAMETER :: routineN = 'pw_integral_ab_fast', 01950 routineP = moduleN//':'//routineN 01951 01952 LOGICAL :: failure 01953 01954 failure = .FALSE. 01955 IF ( pw1%pw_grid %id_nr /= pw2%pw_grid %id_nr ) THEN 01956 CALL stop_program(routineN,moduleN,__LINE__,"Grids incompatible") 01957 END IF 01958 01959 ! since the return value is real, only do accurate sum on the real bit ? 01960 IF ( pw1%in_use == REALDATA3D .AND. pw2%in_use == REALDATA3D ) THEN 01961 integral_value = SUM ( pw1%cr3d ( :, :, : ) & 01962 * pw2%cr3d ( :, :, : ) ) 01963 ELSE IF ( pw1%in_use == REALDATA3D & 01964 .AND. pw2%in_use == COMPLEXDATA3D ) THEN 01965 integral_value = SUM ( pw1%cr3d ( :, :, : ) & 01966 * pw2%cc3d ( :, :, : ) ) 01967 ELSE IF ( pw1%in_use == COMPLEXDATA3D & 01968 .AND. pw2%in_use == REALDATA3D ) THEN 01969 integral_value = SUM ( pw1%cc3d ( :, :, : ) & 01970 * pw2%cr3d ( :, :, : ) ) 01971 ELSE IF ( pw1%in_use == COMPLEXDATA3D & 01972 .AND. pw2%in_use == COMPLEXDATA3D ) THEN 01973 integral_value = SUM ( CONJG ( pw1%cc3d ( :, :, : ) ) & 01974 * pw2%cc3d ( :, :, : ) ) 01975 01976 ELSE IF ( pw1%in_use == REALDATA1D & 01977 .AND. pw2%in_use == REALDATA1D ) THEN 01978 integral_value = SUM ( pw1%cr (:) * pw2%cr (:) ) 01979 ELSE IF ( pw1%in_use == REALDATA1D & 01980 .AND. pw2%in_use == COMPLEXDATA1D ) THEN 01981 integral_value = SUM ( pw1%cr (:) * pw2%cc (:) ) 01982 ELSE IF ( pw1%in_use == COMPLEXDATA1D & 01983 .AND. pw2%in_use == REALDATA1D ) THEN 01984 integral_value = SUM ( pw1%cc (:) * pw2%cr (:) ) 01985 ELSE IF ( pw1%in_use == COMPLEXDATA1D & 01986 .AND. pw2%in_use == COMPLEXDATA1D ) THEN 01987 integral_value = SUM ( CONJG ( pw1%cc (:) ) * pw2%cc (:) ) 01988 ELSE 01989 CALL stop_program(routineN,moduleN,__LINE__,"No possible DATA") 01990 END IF 01991 01992 IF ( pw1%in_use == REALDATA3D .OR. pw1%in_use == COMPLEXDATA3D ) THEN 01993 integral_value = integral_value * pw1%pw_grid%dvol 01994 ELSE 01995 integral_value = integral_value * pw1%pw_grid%vol 01996 ENDIF 01997 IF ( pw1%in_use == COMPLEXDATA1D ) THEN 01998 IF ( pw1%pw_grid%grid_span == HALFSPACE ) THEN 01999 integral_value = 2.0_dp * integral_value 02000 IF ( pw1%pw_grid%have_g0 ) integral_value = integral_value - & 02001 CONJG ( pw1%cc ( 1 ) ) * pw2%cc ( 1 ) 02002 END IF 02003 END IF 02004 02005 IF ( pw1%pw_grid%para%mode == PW_MODE_DISTRIBUTED ) & 02006 CALL mp_sum ( integral_value, pw1%pw_grid%para%group ) 02007 02008 END FUNCTION pw_integral_ab_fast 02009 02010 ! ***************************************************************************** 02011 FUNCTION pw_integral_aa ( pw1, flag, error) RESULT ( integral_value ) 02012 02013 TYPE(pw_type), INTENT(IN) :: pw1 02014 INTEGER, INTENT(IN), OPTIONAL :: flag 02015 TYPE(cp_error_type), INTENT(inout) :: error 02016 REAL(KIND=dp) :: integral_value 02017 02018 CHARACTER(len=*), PARAMETER :: routineN = 'pw_integral_aa', 02019 routineP = moduleN//':'//routineN 02020 02021 LOGICAL :: failure 02022 02023 failure = .FALSE. 02024 IF ( PRESENT ( flag ) ) THEN 02025 IF ( flag == SQUARE ) THEN 02026 IF ( pw1%in_use == REALDATA3D ) THEN 02027 integral_value = accurate_sum ( pw1%cr3d ( :, :, : )** 2 ) 02028 ELSE IF ( pw1%in_use == COMPLEXDATA3D ) THEN 02029 integral_value = accurate_sum ( CONJG ( pw1%cc3d ( :, :, : ) ) & 02030 * pw1%cc3d ( :, :, : ) ) 02031 ELSE IF ( pw1%in_use == REALDATA1D ) THEN 02032 integral_value = accurate_sum ( pw1%cr (:) ** 2 ) 02033 ELSE IF ( pw1%in_use == COMPLEXDATA1D ) THEN 02034 integral_value = accurate_sum ( CONJG ( pw1%cc (:) ) & 02035 * pw1%cc (:) ) 02036 IF ( pw1%pw_grid%grid_span == HALFSPACE ) THEN 02037 integral_value = 2.0_dp * integral_value 02038 IF ( pw1%pw_grid%have_g0 ) integral_value = integral_value - & 02039 CONJG ( pw1%cc ( 1 ) ) * pw1%cc ( 1 ) 02040 END IF 02041 ELSE 02042 CALL stop_program(routineN,moduleN,__LINE__,"No possible SQUARE DATA") 02043 END IF 02044 02045 ELSE IF ( flag == SQUAREROOT ) THEN 02046 CALL stop_program(routineN,moduleN,__LINE__,"No SQUAREROOT defined") 02047 END IF 02048 02049 ELSE 02050 IF ( pw1%in_use == REALDATA3D ) THEN 02051 integral_value = accurate_sum ( pw1%cr3d ( :, :, : ) ) 02052 ELSE IF ( pw1%in_use == COMPLEXDATA3D ) THEN 02053 integral_value = accurate_sum ( pw1%cc3d ( :, :, : ) ) 02054 ELSE IF ( pw1%in_use == REALDATA1D ) THEN 02055 integral_value = accurate_sum ( pw1%cr (:) ) 02056 ELSE IF ( pw1%in_use == COMPLEXDATA1D ) THEN 02057 integral_value = accurate_sum ( pw1%cc (:) ) 02058 ELSE 02059 CALL stop_program(routineN,moduleN,__LINE__,"No possible DATA") 02060 END IF 02061 END IF 02062 02063 IF ( pw1%in_use == REALDATA3D .OR. pw1%in_use == COMPLEXDATA3D ) THEN 02064 integral_value = integral_value * pw1%pw_grid%dvol 02065 ELSE 02066 integral_value = integral_value * pw1%pw_grid%vol 02067 END IF 02068 02069 IF ( pw1%pw_grid%para%mode == PW_MODE_DISTRIBUTED ) & 02070 CALL mp_sum ( integral_value, pw1%pw_grid%para%group ) 02071 02072 END FUNCTION pw_integral_aa 02073 02074 ! ***************************************************************************** 02075 FUNCTION pw_integral_aa_fast ( pw1, flag, error) RESULT ( integral_value ) 02076 02077 TYPE(pw_type), INTENT(IN) :: pw1 02078 INTEGER, INTENT(IN), OPTIONAL :: flag 02079 TYPE(cp_error_type), INTENT(inout) :: error 02080 REAL(KIND=dp) :: integral_value 02081 02082 CHARACTER(len=*), PARAMETER :: routineN = 'pw_integral_aa_fast', 02083 routineP = moduleN//':'//routineN 02084 02085 LOGICAL :: failure 02086 02087 failure = .FALSE. 02088 IF ( PRESENT ( flag ) ) THEN 02089 IF ( flag == SQUARE ) THEN 02090 IF ( pw1%in_use == REALDATA3D ) THEN 02091 integral_value = SUM ( pw1%cr3d ( :, :, : )** 2 ) 02092 ELSE IF ( pw1%in_use == COMPLEXDATA3D ) THEN 02093 integral_value = SUM ( CONJG ( pw1%cc3d ( :, :, : ) ) & 02094 * pw1%cc3d ( :, :, : ) ) 02095 ELSE IF ( pw1%in_use == REALDATA1D ) THEN 02096 integral_value = SUM ( pw1%cr (:) ** 2 ) 02097 ELSE IF ( pw1%in_use == COMPLEXDATA1D ) THEN 02098 integral_value = SUM ( CONJG ( pw1%cc (:) ) & 02099 * pw1%cc (:) ) 02100 IF ( pw1%pw_grid%grid_span == HALFSPACE ) THEN 02101 integral_value = 2.0_dp * integral_value 02102 IF ( pw1%pw_grid%have_g0 ) integral_value = integral_value - & 02103 CONJG ( pw1%cc ( 1 ) ) * pw1%cc ( 1 ) 02104 END IF 02105 ELSE 02106 CALL stop_program(routineN,moduleN,__LINE__,"No possible SQUARE DATA") 02107 END IF 02108 02109 ELSE IF ( flag == SQUAREROOT ) THEN 02110 CALL stop_program(routineN,moduleN,__LINE__,"No SQUAREROOT defined") 02111 END IF 02112 02113 ELSE 02114 IF ( pw1%in_use == REALDATA3D ) THEN 02115 integral_value = SUM ( pw1%cr3d ( :, :, : ) ) 02116 ELSE IF ( pw1%in_use == COMPLEXDATA3D ) THEN 02117 integral_value = SUM ( pw1%cc3d ( :, :, : ) ) 02118 ELSE IF ( pw1%in_use == REALDATA1D ) THEN 02119 integral_value = SUM ( pw1%cr (:) ) 02120 ELSE IF ( pw1%in_use == COMPLEXDATA1D ) THEN 02121 integral_value = SUM ( pw1%cc (:) ) 02122 ELSE 02123 CALL stop_program(routineN,moduleN,__LINE__,"No possible DATA") 02124 END IF 02125 END IF 02126 02127 IF ( pw1%in_use == REALDATA3D .OR. pw1%in_use == COMPLEXDATA3D ) THEN 02128 integral_value = integral_value * pw1%pw_grid%dvol 02129 ELSE 02130 integral_value = integral_value * pw1%pw_grid%vol 02131 END IF 02132 02133 IF ( pw1%pw_grid%para%mode == PW_MODE_DISTRIBUTED ) & 02134 CALL mp_sum ( integral_value, pw1%pw_grid%para%group ) 02135 02136 END FUNCTION pw_integral_aa_fast 02137 02138 ! ***************************************************************************** 02139 FUNCTION pw_integral_a2b ( pw1, pw2, error) RESULT ( integral_value ) 02140 02141 TYPE(pw_type), INTENT(IN) :: pw1, pw2 02142 TYPE(cp_error_type), INTENT(inout) :: error 02143 REAL(KIND=dp) :: integral_value 02144 02145 CHARACTER(len=*), PARAMETER :: routineN = 'pw_integral_a2b', 02146 routineP = moduleN//':'//routineN 02147 02148 LOGICAL :: failure 02149 02150 failure = .FALSE. 02151 IF ( pw1%pw_grid %id_nr /= pw2%pw_grid %id_nr ) THEN 02152 CALL stop_program(routineN,moduleN,__LINE__,"Grids incompatible") 02153 END IF 02154 IF ( pw1%in_use == REALDATA1D .AND. & 02155 pw2%in_use == REALDATA1D ) THEN 02156 integral_value = accurate_sum ( pw1%cr (:) * pw2%cr (:) & 02157 * pw1%pw_grid%gsq (:) ) 02158 ELSE IF ( pw1%in_use == COMPLEXDATA1D .AND. & 02159 pw2%in_use == COMPLEXDATA1D ) THEN 02160 integral_value = accurate_sum ( REAL ( CONJG ( pw1%cc (:) ) & * pw2%cc (:) ,KIND=dp) * pw1%pw_grid%gsq (:)) 02161 IF ( pw1%pw_grid%grid_span == HALFSPACE ) THEN 02162 integral_value = 2.0_dp * integral_value 02163 END IF 02164 ELSE 02165 CALL stop_program(routineN,moduleN,__LINE__,"No possible DATA") 02166 END IF 02167 02168 IF ( pw1%in_use == REALDATA3D .OR. pw1%in_use == COMPLEXDATA3D ) THEN 02169 integral_value = integral_value * pw1%pw_grid%dvol 02170 ELSE 02171 integral_value = integral_value * pw1%pw_grid%vol 02172 END IF 02173 02174 IF ( pw1%pw_grid%para%mode == PW_MODE_DISTRIBUTED ) & 02175 CALL mp_sum ( integral_value, pw1%pw_grid%para%group ) 02176 02177 END FUNCTION pw_integral_a2b 02178 02179 ! ***************************************************************************** 02187 SUBROUTINE pw_structure_factor ( sf, r, error) 02188 02189 TYPE(pw_type), INTENT(INOUT) :: sf 02190 REAL(KIND=dp), DIMENSION(:), INTENT(in) :: r 02191 TYPE(cp_error_type), INTENT(inout) :: error 02192 02193 CHARACTER(len=*), PARAMETER :: routineN = 'pw_structure_factor', 02194 routineP = moduleN//':'//routineN 02195 02196 INTEGER :: cnt, handle, ig 02197 LOGICAL :: failure 02198 REAL(KIND=dp) :: arg, flop 02199 02200 failure = .FALSE. 02201 CALL timeset(routineN,handle) 02202 CPPrecondition(sf%ref_count>0,cp_failure_level,routineP,error,failure) 02203 flop = 0.0_dp 02204 02205 IF ( sf%in_space == RECIPROCALSPACE .AND. & 02206 sf%in_use == COMPLEXDATA1D ) THEN 02207 02208 cnt = SIZE ( sf%cc ) 02209 02210 !$omp parallel do private (ig,arg) default(none) shared(sf,r,cnt) 02211 DO ig = 1, cnt 02212 arg = DOT_PRODUCT(sf%pw_grid%g(:,ig),r) 02213 sf%cc(ig)=CMPLX(COS(arg),-SIN(arg),KIND=dp) 02214 END DO 02215 flop = flop + 7 * cnt 02216 ELSE 02217 02218 CALL stop_program(routineN,moduleN,__LINE__,"No suitable data field") 02219 02220 END IF 02221 02222 flop = flop * 1.e-6_dp 02223 CALL timestop(handle) 02224 02225 END SUBROUTINE pw_structure_factor 02226 02227 ! ***************************************************************************** 02228 FUNCTION pw_integrate_function(fun,isign,oprt,error) RESULT(total_fun) 02229 02230 TYPE(pw_type), INTENT(IN) :: fun 02231 INTEGER, INTENT(IN), OPTIONAL :: isign 02232 CHARACTER(len=*), INTENT(IN), OPTIONAL :: oprt 02233 TYPE(cp_error_type), INTENT(inout) :: error 02234 REAL(KIND=dp) :: total_fun 02235 02236 CHARACTER(len=*), PARAMETER :: routineN = 'pw_integrate_function', 02237 routineP = moduleN//':'//routineN 02238 02239 INTEGER :: iop 02240 LOGICAL :: failure 02241 02242 failure = .FALSE. 02243 iop = 0 02244 IF ( PRESENT(oprt) ) THEN 02245 SELECT CASE (oprt) 02246 CASE ("ABS","abs") 02247 iop = 1 02248 CASE DEFAULT 02249 CALL stop_program(routineN,moduleN,__LINE__,"Unknown operator") 02250 END SELECT 02251 END IF 02252 02253 total_fun = 0._dp 02254 02255 IF (fun%in_space == REALSPACE) THEN 02256 IF (fun%in_use == REALDATA3D) THEN 02257 ! do reduction using maximum accuracy 02258 IF (iop==1) THEN 02259 total_fun = fun%pw_grid%dvol*accurate_sum(ABS(fun%cr3d)) 02260 ELSE 02261 total_fun = fun%pw_grid%dvol*accurate_sum(fun%cr3d) 02262 END IF 02263 ELSE 02264 CALL stop_program(routineN,moduleN,__LINE__,& 02265 "In_space/in_use combination not implemented") 02266 END IF 02267 ELSEIF (fun%in_space == RECIPROCALSPACE) THEN 02268 IF (iop==1) & 02269 CALL stop_program(routineN,moduleN,__LINE__,& 02270 "Operator ABS not implemented") 02271 IF (fun%in_use == COMPLEXDATA3D) THEN 02272 total_fun = fun%pw_grid%vol*fun%cc3d( & 02273 fun%pw_grid%bounds(1,1), & 02274 fun%pw_grid%bounds(1,2), & 02275 fun%pw_grid%bounds(1,3)) 02276 ELSEIF (fun%in_use == COMPLEXDATA1D) THEN 02277 IF ( fun%pw_grid%have_g0 ) total_fun = fun%pw_grid%vol*fun%cc(1) 02278 ELSE 02279 CALL stop_program(routineN,moduleN,__LINE__,& 02280 "In_space/in_use combination not implemented") 02281 END IF 02282 ELSE 02283 CALL stop_program(routineN,moduleN,__LINE__,"No space defined") 02284 END IF 02285 IF (fun%pw_grid%para%mode /= PW_MODE_LOCAL) THEN 02286 CALL mp_sum(total_fun,fun%pw_grid%para%group) 02287 END IF 02288 IF ( PRESENT(isign) ) THEN 02289 total_fun = total_fun * SIGN(1._dp,REAL(isign,dp)) 02290 END IF 02291 02292 END FUNCTION pw_integrate_function 02293 02294 END MODULE pw_methods 02295 02296
1.7.3