|
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 ! ***************************************************************************** 00020 MODULE pw_poisson_types 00021 USE bessel_lib, ONLY: bessj0,& 00022 bessj1,& 00023 bessk0,& 00024 bessk1 00025 USE cell_types, ONLY: cell_release,& 00026 cell_type 00027 USE f77_blas 00028 USE input_constants, ONLY: & 00029 ANALYTIC0D, ANALYTIC1D, ANALYTIC2D, MT0D, MT1D, MT2D, MULTIPOLE0D, & 00030 PERIODIC3D, do_ewald_none, do_ewald_spme, use_analytic, use_mt, & 00031 use_multipole, use_none, use_perd_none, use_perd_x, use_perd_xy, & 00032 use_perd_xyz, use_perd_xz, use_perd_y, use_perd_yz, use_perd_z, & 00033 use_periodic 00034 USE input_section_types, ONLY: section_vals_get_subs_vals,& 00035 section_vals_release,& 00036 section_vals_type,& 00037 section_vals_val_get 00038 USE kinds, ONLY: dp 00039 USE mathconstants, ONLY: fourpi,& 00040 twopi 00041 USE mt_util, ONLY: MTin_create_screen_fn 00042 USE ps_wavelet_types, ONLY: ps_wavelet_release,& 00043 ps_wavelet_type 00044 USE pw_grid_types, ONLY: pw_grid_type 00045 USE pw_grids, ONLY: pw_grid_release 00046 USE pw_pool_types, ONLY: pw_pool_create_pw,& 00047 pw_pool_give_back_pw,& 00048 pw_pool_p_type,& 00049 pw_pool_type,& 00050 pw_pools_dealloc 00051 USE pw_types, ONLY: COMPLEXDATA1D,& 00052 REALDATA1D,& 00053 RECIPROCALSPACE,& 00054 pw_release,& 00055 pw_type 00056 #include "cp_common_uses.h" 00057 00058 IMPLICIT NONE 00059 PRIVATE 00060 00061 LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE. 00062 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_poisson_types' 00063 00064 PUBLIC :: pw_poisson_type 00065 PUBLIC :: pw_poisson_create, pw_poisson_retain, & 00066 pw_poisson_release 00067 PUBLIC :: greens_fn_type,pw_green_create,& 00068 pw_green_retain,& 00069 pw_green_release 00070 00071 INTEGER, SAVE, PRIVATE :: last_greens_fn_id_nr=0 00072 INTEGER, SAVE, PRIVATE :: last_poisson_id=0 00073 00074 ! ***************************************************************************** 00078 TYPE pw_poisson_type 00079 INTEGER :: ref_count, id_nr 00080 INTEGER :: pw_level 00081 INTEGER :: method 00082 INTEGER :: used_grid 00083 LOGICAL :: rebuild 00084 TYPE ( greens_fn_type ), POINTER :: green_fft 00085 TYPE (ps_wavelet_type ),POINTER :: wavelet 00086 TYPE (section_vals_type), POINTER :: parameters 00087 TYPE (cell_type), POINTER :: cell 00088 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools 00089 TYPE ( pw_grid_type), POINTER :: mt_super_ref_pw_grid 00090 END TYPE pw_poisson_type 00091 00092 ! ***************************************************************************** 00096 TYPE greens_fn_type 00097 INTEGER :: method 00098 INTEGER :: special_dimension 00099 INTEGER :: id_nr 00100 INTEGER :: ref_count 00101 REAL (KIND=dp) :: radius 00102 REAL (KIND=dp) :: MT_alpha 00103 REAL (KIND=dp) :: MT_rel_cutoff 00104 REAL (KIND=dp) :: slab_size 00105 REAL (KIND=dp) :: alpha 00106 LOGICAL :: p3m 00107 INTEGER :: p3m_order 00108 REAL (KIND=dp) :: p3m_alpha 00109 REAL (KIND=dp), DIMENSION ( :, : ), POINTER :: p3m_coeff 00110 REAL (KIND=dp), DIMENSION ( :, : ), POINTER :: p3m_bm2 00111 LOGICAL :: sr_screening 00112 REAL (KIND=dp) :: sr_alpha 00113 REAL (KIND=dp) :: sr_rc 00114 TYPE ( pw_type ), POINTER :: influence_fn 00115 TYPE ( pw_type ), POINTER :: screen_fn 00116 TYPE ( pw_type ), POINTER :: p3m_charge 00117 END TYPE greens_fn_type 00118 00119 CONTAINS 00120 00121 ! ***************************************************************************** 00126 SUBROUTINE pw_green_create ( green, poisson_section, cell, pw_pool, & 00127 mt_super_ref_pw_grid, error ) 00128 TYPE(greens_fn_type), POINTER :: green 00129 TYPE(section_vals_type), POINTER :: poisson_section 00130 TYPE(cell_type), POINTER :: cell 00131 TYPE(pw_pool_type), POINTER :: pw_pool 00132 TYPE(pw_grid_type), POINTER :: mt_super_ref_pw_grid 00133 TYPE(cp_error_type), INTENT(inout) :: error 00134 00135 CHARACTER(len=*), PARAMETER :: routineN = 'pw_green_create', 00136 routineP = moduleN//':'//routineN 00137 INTEGER, PARAMETER :: nr = 250 00138 00139 INTEGER :: dim, ewald_kind, i, ig, iz, 00140 my_per, my_val, n, nz, stat 00141 INTEGER, DIMENSION(3) :: perd 00142 LOGICAL :: failure 00143 REAL(KIND=dp) :: g2, g3d, gg, gxy, j0g, j1g, 00144 k0g, k1g, rlength, zlength 00145 REAL(KIND=dp), DIMENSION(3) :: abc 00146 TYPE(pw_grid_type), POINTER :: grid 00147 TYPE(pw_type), POINTER :: gf 00148 TYPE(section_vals_type), POINTER :: ewald_section, mt_section 00149 00150 failure = .FALSE. 00151 NULLIFY(ewald_section) 00152 CPPrecondition(ASSOCIATED(cell),cp_failure_level,routineP,error,failure) 00153 CPPrecondition(.NOT.(ASSOCIATED(green)),cp_failure_level,routineP,error,failure) 00154 IF (.NOT.failure) THEN 00155 ALLOCATE(green, stat=stat) 00156 CPPostcondition(stat == 0,cp_fatal_level,routineP,error,failure) 00157 END IF 00158 IF (.NOT.failure) THEN 00159 green%p3m=.FALSE. 00160 green%special_dimension = 0 00161 green%radius = 0.0_dp 00162 green%slab_size = 0.0_dp 00163 green%alpha = 0.0_dp 00164 green%method=PERIODIC3D 00165 last_greens_fn_id_nr = last_greens_fn_id_nr+1 00166 green % id_nr = last_greens_fn_id_nr 00167 green%ref_count=1 00168 green%MT_alpha=1.0_dp 00169 green%MT_rel_cutoff=1.0_dp 00170 green%p3m=.FALSE. 00171 green%p3m_order=0 00172 green%p3m_alpha=0.0_dp 00173 green%sr_screening=.FALSE. 00174 green%sr_alpha=1.0_dp 00175 green%sr_rc=0.0_dp 00176 00177 NULLIFY (green%influence_fn,green%p3m_charge) 00178 NULLIFY (green%p3m_coeff,green%p3m_bm2) 00179 NULLIFY (green%screen_fn) 00180 00181 CALL section_vals_val_get(poisson_section,"PERIODIC",i_val=my_per,error=error) 00182 ewald_section => section_vals_get_subs_vals(poisson_section,& 00183 "EWALD", can_return_null=.TRUE., error=error) 00184 ewald_kind=do_ewald_none 00185 IF ( ASSOCIATED(ewald_section) ) THEN 00186 CALL section_vals_val_get(ewald_section,"EWALD_TYPE",i_val=ewald_kind,& 00187 error=error) 00188 END IF 00189 00190 !CPPrecondition(cell%orthorhombic,cp_failure_level,routineP,error,failure) 00191 DO i=1,3 00192 abc(i)=cell%hmat(i,i) 00193 END DO 00194 SELECT CASE(my_per) 00195 CASE(use_perd_x) 00196 perd = (/1,0,0/) 00197 CASE(use_perd_y) 00198 perd = (/0,1,0/) 00199 CASE(use_perd_z) 00200 perd = (/0,0,1/) 00201 CASE(use_perd_xy) 00202 perd = (/1,1,0/) 00203 CASE(use_perd_xz) 00204 perd = (/1,0,1/) 00205 CASE(use_perd_yz) 00206 perd = (/0,1,1/) 00207 CASE(use_perd_xyz) 00208 perd = (/1,1,1/) 00209 CASE(use_perd_none) 00210 perd = (/0,0,0/) 00211 CASE DEFAULT 00212 CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) 00213 END SELECT 00214 ! check for consistent use of periodicity (cell <-> Poisson solver) 00215 !CPPostcondition(ALL(perd == cell%perd),cp_fatal_level,routineP,error,failure) 00216 00217 dim = COUNT(perd == 1) 00218 CALL section_vals_val_get(poisson_section,"POISSON_SOLVER",i_val=my_val,error=error) 00219 SELECT CASE (my_val) 00220 CASE (use_periodic) 00221 green%method = PERIODIC3D 00222 IF (dim /= 3) THEN 00223 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 00224 "Illegal combination of periodicity and Poisson solver periodic3d",& 00225 error=error,failure=failure) 00226 END IF 00227 CASE (use_multipole) 00228 green%method = MULTIPOLE0D 00229 IF (dim /= 0) THEN 00230 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 00231 "Illegal combination of periodicity and Poisson solver mulipole0d",& 00232 error=error,failure=failure) 00233 END IF 00234 CASE (use_analytic) 00235 SELECT CASE (dim) 00236 CASE (0) 00237 green%method = ANALYTIC0D 00238 green%radius = 0.5_dp*MINVAL(abc) 00239 CASE (1) 00240 green%method = ANALYTIC1D 00241 green%special_dimension = MAXLOC(perd(1:3),1) 00242 green%radius = MAXVAL(abc(1:3)) 00243 DO i=1,3 00244 IF (i == green%special_dimension) CYCLE 00245 green%radius = MIN(green%radius,0.5_dp*abc(i)) 00246 END DO 00247 CASE (2) 00248 green%method = ANALYTIC2D 00249 i = MINLOC(perd,1) 00250 green%special_dimension = i 00251 green%slab_size = abc(i) 00252 CASE (3) 00253 green%method = PERIODIC3D 00254 CASE DEFAULT 00255 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00256 END SELECT 00257 CASE (use_mt) 00258 mt_section => section_vals_get_subs_vals(poisson_section,"MT",error=error) 00259 CALL section_vals_val_get(mt_section,"REL_CUTOFF",r_val=green%MT_rel_cutoff,error=error) 00260 CALL section_vals_val_get(mt_section,"ALPHA",r_val=green%MT_alpha,error=error) 00261 green%MT_alpha=green%MT_alpha/MINVAL(abc) 00262 SELECT CASE (dim) 00263 CASE (0) 00264 green%method = MT0D 00265 green%radius = 0.5_dp*MINVAL(abc) 00266 CASE (1) 00267 green%method = MT1D 00268 green%special_dimension = MAXLOC(perd(1:3),1) 00269 green%radius = MAXVAL(abc(1:3)) 00270 DO i=1,3 00271 IF (i == green%special_dimension) CYCLE 00272 green%radius = MIN(green%radius,0.5_dp*abc(i)) 00273 END DO 00274 CASE (2) 00275 green%method = MT2D 00276 i = MINLOC(perd,1) 00277 green%special_dimension = i 00278 green%slab_size = abc(i) 00279 CASE (3) 00280 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 00281 "Illegal combination of periodicity and Poisson solver (MT)",& 00282 error=error,failure=failure) 00283 CASE DEFAULT 00284 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00285 END SELECT 00286 CASE DEFAULT 00287 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 00288 "An unknown Poisson solver was specified",error,failure) 00289 END SELECT 00290 ! allocate influence function,... 00291 SELECT CASE ( green % method ) 00292 CASE ( PERIODIC3D, ANALYTIC2D, ANALYTIC1D, ANALYTIC0D, MT2D, MT1D, MT0D, MULTIPOLE0D) 00293 CALL pw_pool_create_pw ( pw_pool, green % influence_fn,& 00294 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE ,error=error) 00295 IF (ewald_kind==do_ewald_spme) THEN 00296 green % p3m = .TRUE. 00297 CALL section_vals_val_get(ewald_section,"o_spline",& 00298 i_val=green % p3m_order,error=error) 00299 CALL section_vals_val_get(ewald_section,"alpha",& 00300 r_val=green%p3m_alpha,error=error) 00301 n=green % p3m_order 00302 ALLOCATE ( green%p3m_coeff ( -(n-1):n-1, 0:n-1 ), stat = stat ) 00303 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 00304 CALL spme_coeff_calculate ( n, green%p3m_coeff ) 00305 CALL pw_pool_create_pw ( pw_pool, green % p3m_charge, use_data=REALDATA1D, & 00306 in_space=RECIPROCALSPACE,error=error) 00307 CALL influence_factor ( green ,error=error) 00308 CALL calc_p3m_charge(green) 00309 ELSE 00310 green % p3m = .FALSE. 00311 END IF 00312 ! 00313 SELECT CASE(green%method) 00314 CASE(MT0D,MT1D,MT2D) 00315 CALL MTin_create_screen_fn(green%screen_fn,pw_pool=pw_pool,method=green%method,& 00316 alpha=green%MT_alpha, & 00317 special_dimension=green%special_dimension, slab_size=green%slab_size, & 00318 super_ref_pw_grid=mt_super_ref_pw_grid, error=error) 00319 END SELECT 00320 CASE DEFAULT 00321 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00322 END SELECT 00323 00324 ! initialize influence function 00325 gf => green % influence_fn 00326 grid => green % influence_fn % pw_grid 00327 SELECT CASE ( green%method ) 00328 CASE ( PERIODIC3D, MULTIPOLE0D ) 00329 DO ig = grid % first_gne0, grid % ngpts_cut_local 00330 g2 = grid % gsq ( ig ) 00331 gf % cc ( ig ) = fourpi / g2 00332 END DO 00333 IF ( grid % have_g0 ) gf % cc ( 1 ) = 0.0_dp 00334 00335 CASE ( ANALYTIC2D ) 00336 00337 iz = green % special_dimension ! iz is the direction with NO PBC 00338 zlength = green % slab_size ! zlength is the thickness of the cell 00339 DO ig = grid % first_gne0, grid % ngpts_cut_local 00340 nz = grid % g_hat ( iz, ig ) 00341 g2 = grid % gsq ( ig ) 00342 g3d = fourpi / g2 00343 gg = 0.5_dp * SQRT ( g2 ) 00344 gf % cc ( ig ) = g3d * ( 1.0_dp - (-1.0_dp)**nz * EXP ( - gg * zlength ) ) 00345 END DO 00346 IF ( grid % have_g0 ) gf % cc ( 1 ) = 0.0_dp 00347 00348 CASE ( ANALYTIC1D ) 00349 00350 ! iz is the direction of the PBC ( can be 1,2,3 -> x,y,z ) 00351 iz = green % special_dimension 00352 ! rlength is the radius of the tube 00353 rlength = green % radius 00354 DO ig = grid % first_gne0, grid % ngpts_cut_local 00355 g2 = grid % gsq ( ig ) 00356 g3d = fourpi / g2 00357 gxy = SQRT ( g2 - grid % g(iz,ig) * grid % g(iz,ig) ) 00358 j0g = bessj0 ( rlength * gxy ) 00359 j1g = bessj1 ( rlength * gxy ) 00360 k0g = bessk0 ( rlength * grid % g(iz,ig) ) 00361 k1g = bessk1 ( rlength * grid % g(iz,ig) ) 00362 gf % cc ( ig ) = g3d * ( 1.0_dp - rlength * & 00363 ( gxy * j1g * k0g - grid % g(iz,ig) * j0g * k1g ) ) 00364 END DO 00365 IF ( grid % have_g0 ) gf % cc ( 1 ) = 0.0_dp 00366 00367 CASE ( ANALYTIC0D ) 00368 00369 rlength = green % radius ! rlength is the radius of the sphere 00370 DO ig = grid % first_gne0, grid % ngpts_cut_local 00371 g2 = grid % gsq ( ig ) 00372 gg = SQRT ( g2 ) 00373 g3d = fourpi / g2 00374 gf % cc ( ig ) = g3d * ( 1.0_dp - COS ( rlength * gg ) ) 00375 END DO 00376 IF ( grid % have_g0 ) & 00377 gf % cc ( 1 ) = 0.5_dp * fourpi * rlength * rlength 00378 00379 CASE ( MT2D, MT1D, MT0D ) 00380 DO ig = grid % first_gne0, grid % ngpts_cut_local 00381 g2 = grid % gsq ( ig ) 00382 g3d = fourpi / g2 00383 gf%cc ( ig ) = g3d + green%screen_fn%cc(ig) 00384 END DO 00385 IF ( grid % have_g0 ) & 00386 gf%cc(1) = green%screen_fn%cc(1) 00387 CASE DEFAULT 00388 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00389 END SELECT 00390 00391 END IF 00392 END SUBROUTINE pw_green_create 00393 00394 ! ***************************************************************************** 00398 SUBROUTINE pw_green_retain(gftype,error) 00399 TYPE(greens_fn_type), POINTER :: gftype 00400 TYPE(cp_error_type), INTENT(inout) :: error 00401 00402 CHARACTER(len=*), PARAMETER :: routineN = 'pw_green_retain', 00403 routineP = moduleN//':'//routineN 00404 00405 LOGICAL :: failure 00406 00407 failure=.FALSE. 00408 00409 CPPrecondition(ASSOCIATED(gftype),cp_failure_level,routineP,error,failure) 00410 IF (.NOT. failure) THEN 00411 CPPrecondition(gftype%ref_count>0,cp_failure_level,routineP,error,failure) 00412 gftype%ref_count=gftype%ref_count+1 00413 END IF 00414 END SUBROUTINE pw_green_retain 00415 00416 ! ***************************************************************************** 00423 SUBROUTINE pw_green_release ( gftype, pw_pool, error ) 00424 TYPE(greens_fn_type), POINTER :: gftype 00425 TYPE(pw_pool_type), OPTIONAL, POINTER :: pw_pool 00426 TYPE(cp_error_type), INTENT(inout) :: error 00427 00428 CHARACTER(len=*), PARAMETER :: routineN = 'pw_green_release', 00429 routineP = moduleN//':'//routineN 00430 00431 INTEGER :: stat 00432 LOGICAL :: can_give_back, failure 00433 00434 failure = .FALSE. 00435 IF (.NOT.failure) THEN 00436 IF (ASSOCIATED(gftype)) THEN 00437 CPPrecondition(gftype%ref_count>0,cp_failure_level,routineP,error,failure) 00438 gftype%ref_count=gftype%ref_count-1 00439 IF (gftype%ref_count==0) THEN 00440 can_give_back=PRESENT(pw_pool) 00441 IF (can_give_back) can_give_back=ASSOCIATED(pw_pool) 00442 IF (can_give_back) THEN 00443 CALL pw_pool_give_back_pw(pw_pool,gftype%influence_fn,& 00444 accept_non_compatible=.TRUE.,error=error) 00445 CALL pw_pool_give_back_pw(pw_pool,gftype%screen_fn,& 00446 accept_non_compatible=.TRUE.,error=error) 00447 CALL pw_pool_give_back_pw(pw_pool,gftype%p3m_charge,& 00448 accept_non_compatible=.TRUE.,error=error) 00449 ELSE 00450 CALL pw_release(gftype%influence_fn,error=error) 00451 CALL pw_release(gftype%screen_fn,error=error) 00452 CALL pw_release(gftype % p3m_charge, error=error) 00453 END IF 00454 IF (ASSOCIATED(gftype % p3m_bm2)) & 00455 DEALLOCATE ( gftype % p3m_bm2 ) 00456 IF (ASSOCIATED(gftype % p3m_coeff)) & 00457 DEALLOCATE ( gftype % p3m_coeff ) 00458 DEALLOCATE(gftype, stat=stat) 00459 CPPostcondition(stat == 0,cp_failure_level,routineP,error,failure) 00460 END IF 00461 END IF 00462 END IF 00463 NULLIFY(gftype) 00464 END SUBROUTINE pw_green_release 00465 00466 ! ***************************************************************************** 00473 SUBROUTINE influence_factor ( gftype, error ) 00474 TYPE(greens_fn_type), POINTER :: gftype 00475 TYPE(cp_error_type), INTENT(inout) :: error 00476 00477 CHARACTER(len=*), PARAMETER :: routineN = 'influence_factor', 00478 routineP = moduleN//':'//routineN 00479 00480 COMPLEX(KIND=dp) :: b_m, exp_m, sum_m 00481 INTEGER :: dim, ierr, j, k, l, n, pt 00482 INTEGER, DIMENSION(3) :: npts 00483 INTEGER, DIMENSION(:), POINTER :: lb, ub 00484 LOGICAL :: failure 00485 REAL(KIND=dp) :: l_arg, prod_arg, val 00486 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: m_assign 00487 00488 failure=.FALSE. 00489 CPPrecondition(ASSOCIATED(gftype),cp_failure_level,routineP,error,failure) 00490 CPPrecondition(gftype%ref_count>0,cp_failure_level,routineP,error,failure) 00491 n = gftype % p3m_order 00492 00493 ! calculate the assignment function values 00494 00495 lb => gftype % influence_fn % pw_grid % bounds (1, : ) 00496 ub => gftype % influence_fn % pw_grid % bounds (2, : ) 00497 IF (ASSOCIATED(gftype % p3m_bm2)) THEN 00498 IF (LBOUND(gftype % p3m_bm2,2)/=MINVAL(lb(:)).OR.& 00499 UBOUND(gftype % p3m_bm2,2)/=MAXVAL(ub(:))) THEN 00500 DEALLOCATE(gftype % p3m_bm2,stat=ierr) 00501 CPPostcondition(ierr==0,cp_fatal_level,routineP,error,failure) 00502 END IF 00503 END IF 00504 IF (.NOT.ASSOCIATED(gftype % p3m_bm2)) THEN 00505 ALLOCATE ( gftype % p3m_bm2 ( 3, MINVAL(lb(:)):MAXVAL(ub(:)) ), STAT = ierr ) 00506 CPPostcondition(ierr==0,cp_fatal_level,routineP,error,failure) 00507 END IF 00508 00509 ALLOCATE ( m_assign ( 0:n-2 ), STAT = ierr ) 00510 CPPostcondition(ierr==0,cp_fatal_level,routineP,error,failure) 00511 m_assign = 0.0_dp 00512 DO k = 0, n-2 00513 j = -(n-1) + 2 * k 00514 DO l = 0, n-1 00515 l_arg = 0.5_dp ** l 00516 prod_arg = gftype % p3m_coeff ( j, l ) * l_arg 00517 m_assign ( k ) = m_assign ( k ) + prod_arg 00518 END DO 00519 END DO 00520 00521 ! calculate the absolute b values 00522 00523 npts ( : ) = ub ( : ) - lb ( : ) + 1 00524 DO dim = 1, 3 00525 DO pt = lb (dim), ub (dim) 00526 val = twopi * ( REAL ( pt,KIND=dp) / REAL ( npts ( dim ),KIND=dp) ) 00527 exp_m = CMPLX ( COS ( val ), -SIN ( val ),KIND=dp) 00528 sum_m = CMPLX ( 0.0_dp, 0.0_dp,KIND=dp) 00529 DO k = 0, n-2 00530 sum_m = sum_m + m_assign ( k ) * exp_m ** k 00531 END DO 00532 b_m = exp_m ** ( n - 1 ) / sum_m 00533 gftype % p3m_bm2 ( dim, pt ) = SQRT ( REAL ( b_m * CONJG ( b_m ),KIND=dp) ) 00534 END DO 00535 END DO 00536 00537 DEALLOCATE ( m_assign, STAT = ierr ) 00538 CPPostconditionNoFail(ierr==0,cp_warning_level,routineP,error) 00539 END SUBROUTINE influence_factor 00540 00541 ! ***************************************************************************** 00542 SUBROUTINE calc_p3m_charge ( gf ) 00543 00544 TYPE(greens_fn_type), POINTER :: gf 00545 00546 INTEGER :: ig, l, m, n 00547 REAL(KIND=dp) :: arg, novol 00548 REAL(KIND=dp), DIMENSION(:, :), POINTER :: bm2 00549 TYPE(pw_grid_type), POINTER :: grid 00550 TYPE(pw_type), POINTER :: pc 00551 00552 grid => gf % influence_fn % pw_grid 00553 00554 ! check if charge function is consistent with current box volume 00555 00556 pc => gf % p3m_charge 00557 bm2 => gf % p3m_bm2 00558 arg = 1.0_dp / ( 8.0_dp * gf % p3m_alpha ** 2 ) 00559 novol = REAL ( grid % ngpts,KIND=dp) / grid % vol 00560 DO ig = 1, grid % ngpts_cut_local 00561 l = grid % g_hat ( 1, ig ) 00562 m = grid % g_hat ( 2, ig ) 00563 n = grid % g_hat ( 3, ig ) 00564 pc % cr ( ig ) = novol * EXP ( -arg * grid % gsq ( ig ) ) * & 00565 bm2 ( 1, l ) * bm2 ( 2, m ) * bm2 ( 3, n ) 00566 END DO 00567 00568 END SUBROUTINE calc_p3m_charge 00569 00570 ! ***************************************************************************** 00579 SUBROUTINE pw_poisson_create ( poisson_env, error ) 00580 00581 TYPE(pw_poisson_type), POINTER :: poisson_env 00582 TYPE(cp_error_type), INTENT(inout) :: error 00583 00584 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_create', 00585 routineP = moduleN//':'//routineN 00586 00587 INTEGER :: stat 00588 LOGICAL :: failure 00589 00590 failure=.FALSE. 00591 00592 CPPrecondition(.NOT.ASSOCIATED(poisson_env),cp_failure_level,routineP,error,failure) 00593 IF (.NOT.failure) THEN 00594 ALLOCATE(poisson_env,stat=stat) 00595 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 00596 END IF 00597 IF (.NOT.failure) THEN 00598 last_poisson_id=last_poisson_id+1 00599 poisson_env%id_nr=last_poisson_id 00600 poisson_env%ref_count=1 00601 poisson_env%method=use_none 00602 poisson_env%rebuild=.TRUE. 00603 NULLIFY(poisson_env%parameters,poisson_env%green_fft,poisson_env%cell,& 00604 poisson_env%pw_pools,poisson_env%green_fft,& 00605 poisson_env%mt_super_ref_pw_grid,poisson_env%wavelet) 00606 poisson_env%pw_level=-1 00607 poisson_env%ref_count=1 00608 END IF 00609 00610 END SUBROUTINE pw_poisson_create 00611 00612 ! ***************************************************************************** 00618 SUBROUTINE pw_poisson_retain(poisson_env,error) 00619 TYPE(pw_poisson_type), POINTER :: poisson_env 00620 TYPE(cp_error_type), INTENT(inout) :: error 00621 00622 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_retain', 00623 routineP = moduleN//':'//routineN 00624 00625 LOGICAL :: failure 00626 00627 failure=.FALSE. 00628 00629 CPPrecondition(ASSOCIATED(poisson_env),cp_failure_level,routineP,error,failure) 00630 IF (.NOT. failure) THEN 00631 CPPreconditionNoFail(poisson_env%ref_count>0,cp_failure_level,routineP,error) 00632 poisson_env%ref_count=poisson_env%ref_count+1 00633 END IF 00634 END SUBROUTINE pw_poisson_retain 00635 00636 ! ***************************************************************************** 00642 SUBROUTINE pw_poisson_release ( poisson_env, error) 00643 00644 TYPE(pw_poisson_type), POINTER :: poisson_env 00645 TYPE(cp_error_type), INTENT(inout) :: error 00646 00647 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_release', 00648 routineP = moduleN//':'//routineN 00649 00650 INTEGER :: stat 00651 LOGICAL :: failure 00652 00653 failure=.FALSE. 00654 IF (ASSOCIATED(poisson_env)) THEN 00655 CPPrecondition(poisson_env%ref_count>0,cp_failure_level,routineP,error,failure) 00656 poisson_env%ref_count=poisson_env%ref_count-1 00657 IF (poisson_env%ref_count==0) THEN 00658 CALL section_vals_release(poisson_env%parameters,error=error) 00659 IF (ASSOCIATED(poisson_env%pw_pools)) THEN 00660 CALL pw_pools_dealloc(poisson_env%pw_pools,error=error) 00661 END IF 00662 CALL pw_green_release(poisson_env%green_fft,error=error) 00663 CALL cell_release(poisson_env%cell,error=error) 00664 CALL pw_grid_release(poisson_env%mt_super_ref_pw_grid,error=error) 00665 CALL ps_wavelet_release(poisson_env%wavelet,error=error) 00666 DEALLOCATE(poisson_env,stat=stat) 00667 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00668 END IF 00669 END IF 00670 NULLIFY(poisson_env) 00671 00672 END SUBROUTINE pw_poisson_release 00673 00674 ! ***************************************************************************** 00680 SUBROUTINE spme_coeff_calculate ( n, coeff ) 00681 00682 INTEGER, INTENT(IN) :: n 00683 REAL(KIND=dp), 00684 DIMENSION(-(n-1):n-1, 0:n-1), 00685 INTENT(OUT) :: coeff 00686 00687 INTEGER :: i, j, l, m 00688 REAL(KIND=dp) :: b 00689 REAL(KIND=dp), DIMENSION(n, -n:n, 0:n-1) :: a 00690 00691 a = 0.0_dp 00692 a ( 1, 0, 0 ) = 1.0_dp 00693 00694 DO i = 2, n 00695 m = i-1 00696 DO j = -m, m, 2 00697 DO l = 0, m-1 00698 b = ( a ( m, j-1, l ) + & 00699 REAL ( (-1) ** l,KIND=dp) * a ( m, j+1, l ) ) / 00700 REAL ( ( l + 1 ) * 2 ** ( l + 1 ) ,KIND=dp) 00701 a ( i, j, 0 ) = a ( i, j, 0 ) + b 00702 END DO 00703 DO l = 0, m-1 00704 a ( i, j, l+1 ) = ( a ( m, j+1, l ) - & 00705 a ( m, j-1, l ) ) / REAL ( l + 1,KIND=dp) 00706 END DO 00707 END DO 00708 END DO 00709 00710 coeff = 0.0_dp 00711 DO i = 0, n-1 00712 DO j = -(n-1), n-1, 2 00713 coeff ( j, i ) = a ( n, j, i ) 00714 END DO 00715 END DO 00716 00717 END SUBROUTINE spme_coeff_calculate 00718 00719 END MODULE pw_poisson_types
1.7.3