1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief functions related to the poisson solver on regular grids 8!> \par History 9!> greens_fn: JGH (9-Mar-2001) : include influence_fn into 10!> greens_fn_type add cell volume 11!> as indicator for updates 12!> greens_fn: JGH (30-Mar-2001) : Added B-spline routines 13!> pws : JGH (13-Mar-2001) : new pw_poisson_solver, delete 14!> pw_greens_fn 15!> 12.2004 condensed from pws, greens_fn and green_fns, by apsi and JGH, 16!> made thread safe, new input [fawzi] 17!> 14-Mar-2006 : added short range screening function for SE codes 18!> \author fawzi 19! ************************************************************************************************** 20MODULE pw_poisson_types 21 USE bessel_lib, ONLY: bessk0,& 22 bessk1 23 USE dielectric_types, ONLY: dielectric_parameters 24 USE dirichlet_bc_types, ONLY: dirichlet_bc_parameters 25 USE kinds, ONLY: dp 26 USE mathconstants, ONLY: fourpi,& 27 twopi 28 USE mt_util, ONLY: MT0D,& 29 MT1D,& 30 MT2D,& 31 MTin_create_screen_fn 32 USE ps_implicit_types, ONLY: MIXED_BC,& 33 NEUMANN_BC,& 34 ps_implicit_parameters,& 35 ps_implicit_release,& 36 ps_implicit_type 37 USE ps_wavelet_types, ONLY: ps_wavelet_release,& 38 ps_wavelet_type 39 USE pw_grid_types, ONLY: pw_grid_type 40 USE pw_grids, ONLY: pw_grid_release 41 USE pw_pool_types, ONLY: pw_pool_create,& 42 pw_pool_create_pw,& 43 pw_pool_give_back_pw,& 44 pw_pool_p_type,& 45 pw_pool_release,& 46 pw_pool_type,& 47 pw_pools_dealloc 48 USE pw_types, ONLY: COMPLEXDATA1D,& 49 REALDATA1D,& 50 RECIPROCALSPACE,& 51 pw_release,& 52 pw_type 53 USE realspace_grid_types, ONLY: realspace_grid_type,& 54 rs_grid_release 55#include "../base/base_uses.f90" 56 57 IMPLICIT NONE 58 PRIVATE 59 60 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE. 61 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_poisson_types' 62 63 PUBLIC :: pw_poisson_type 64 PUBLIC :: pw_poisson_create, pw_poisson_retain, & 65 pw_poisson_release 66 PUBLIC :: greens_fn_type, pw_green_create, & 67 pw_green_release 68 PUBLIC :: pw_poisson_parameter_type 69 70 INTEGER, SAVE, PRIVATE :: last_greens_fn_id_nr = 0 71 INTEGER, SAVE, PRIVATE :: last_poisson_id = 0 72 73 INTEGER, PARAMETER, PUBLIC :: pw_poisson_none = 0, & 74 pw_poisson_periodic = 1, & 75 pw_poisson_analytic = 2, & 76 pw_poisson_mt = 3, & 77 pw_poisson_hockney = 5, & 78 pw_poisson_multipole = 4, & 79 pw_poisson_wavelet = 6, & 80 pw_poisson_implicit = 7 81 ! EWALD methods 82 INTEGER, PARAMETER, PUBLIC :: do_ewald_none = 1, & 83 do_ewald_ewald = 2, & 84 do_ewald_pme = 3, & 85 do_ewald_spme = 4 86 87 INTEGER, PARAMETER, PUBLIC :: PERIODIC3D = 1000, & 88 ANALYTIC2D = 1001, & 89 ANALYTIC1D = 1002, & 90 ANALYTIC0D = 1003, & 91 HOCKNEY2D = 1201, & 92 HOCKNEY1D = 1202, & 93 HOCKNEY0D = 1203, & 94 MULTIPOLE2D = 1301, & 95 MULTIPOLE1D = 1302, & 96 MULTIPOLE0D = 1303, & 97 PS_IMPLICIT = 1400 98 99! ************************************************************************************************** 100!> \brief parameters for the poisson solver independet of input_section 101!> \author Ole Schuett 102! ************************************************************************************************** 103 TYPE pw_poisson_parameter_type 104 INTEGER :: solver 105 106 INTEGER, DIMENSION(3) :: periodic 107 INTEGER :: ewald_type = do_ewald_none 108 INTEGER :: ewald_o_spline 109 REAL(KIND=dp) :: ewald_alpha 110 111 REAL(KIND=dp) :: mt_rel_cutoff 112 REAL(KIND=dp) :: mt_alpha 113 114 INTEGER :: wavelet_scf_type 115 INTEGER :: wavelet_method 116 INTEGER :: wavelet_special_dimension 117 CHARACTER(LEN=1) :: wavelet_geocode 118 119 LOGICAL :: has_dielectric 120 TYPE(dielectric_parameters) :: dielectric_params 121 TYPE(ps_implicit_parameters) :: ps_implicit_params 122 TYPE(dirichlet_bc_parameters) :: dbc_params 123 END TYPE pw_poisson_parameter_type 124 125! ************************************************************************************************** 126!> \brief environment for the poisson solver 127!> \author fawzi 128! ************************************************************************************************** 129 TYPE pw_poisson_type 130 INTEGER :: ref_count, id_nr 131 INTEGER :: pw_level 132 INTEGER :: method 133 INTEGER :: used_grid 134 LOGICAL :: rebuild 135 TYPE(greens_fn_type), POINTER :: green_fft 136 TYPE(ps_wavelet_type), POINTER :: wavelet 137 TYPE(pw_poisson_parameter_type) :: parameters 138 REAL(KIND=dp), DIMENSION(3, 3) :: cell_hmat = 0.0_dp 139 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools 140 TYPE(pw_grid_type), POINTER :: mt_super_ref_pw_grid 141 TYPE(ps_implicit_type), POINTER :: implicit_env 142 TYPE(pw_grid_type), POINTER :: dct_pw_grid 143 TYPE(realspace_grid_type), POINTER :: diel_rs_grid 144 END TYPE pw_poisson_type 145 146! ************************************************************************************************** 147!> \brief contains all the informations needed by the fft based poisson solvers 148!> \author JGH,Teo,fawzi 149! ************************************************************************************************** 150 TYPE greens_fn_type 151 INTEGER :: method 152 INTEGER :: special_dimension 153 INTEGER :: id_nr 154 INTEGER :: ref_count 155 REAL(KIND=dp) :: radius 156 REAL(KIND=dp) :: MT_alpha 157 REAL(KIND=dp) :: MT_rel_cutoff 158 REAL(KIND=dp) :: slab_size 159 REAL(KIND=dp) :: alpha 160 LOGICAL :: p3m 161 INTEGER :: p3m_order 162 REAL(KIND=dp) :: p3m_alpha 163 REAL(KIND=dp), DIMENSION(:, :), POINTER :: p3m_coeff 164 REAL(KIND=dp), DIMENSION(:, :), POINTER :: p3m_bm2 165 LOGICAL :: sr_screening 166 REAL(KIND=dp) :: sr_alpha 167 REAL(KIND=dp) :: sr_rc 168 TYPE(pw_type), POINTER :: influence_fn 169 TYPE(pw_type), POINTER :: dct_influence_fn 170 TYPE(pw_type), POINTER :: screen_fn 171 TYPE(pw_type), POINTER :: p3m_charge 172 END TYPE greens_fn_type 173 174CONTAINS 175 176! ************************************************************************************************** 177!> \brief Allocates and sets up the green functions for the fft based poisson 178!> solvers 179!> \param green ... 180!> \param poisson_params ... 181!> \param cell_hmat ... 182!> \param pw_pool ... 183!> \param mt_super_ref_pw_grid ... 184!> \param dct_pw_grid ... 185!> \author Fawzi, based on previous functions by JGH and Teo 186! ************************************************************************************************** 187 SUBROUTINE pw_green_create(green, poisson_params, cell_hmat, pw_pool, & 188 mt_super_ref_pw_grid, dct_pw_grid) 189 TYPE(greens_fn_type), POINTER :: green 190 TYPE(pw_poisson_parameter_type), INTENT(IN) :: poisson_params 191 REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: cell_hmat 192 TYPE(pw_pool_type), POINTER :: pw_pool 193 TYPE(pw_grid_type), POINTER :: mt_super_ref_pw_grid, dct_pw_grid 194 195 CHARACTER(len=*), PARAMETER :: routineN = 'pw_green_create', & 196 routineP = moduleN//':'//routineN 197 198 INTEGER :: dim, i, ig, iz, n, nz 199 REAL(KIND=dp) :: g2, g3d, gg, gxy, gz, j0g, j1g, k0g, & 200 k1g, rlength, zlength 201 REAL(KIND=dp), DIMENSION(3) :: abc 202 TYPE(pw_grid_type), POINTER :: dct_grid, grid 203 TYPE(pw_pool_type), POINTER :: pw_pool_xpndd 204 TYPE(pw_type), POINTER :: dct_gf, gf 205 206 CPASSERT(.NOT. (ASSOCIATED(green))) 207 ALLOCATE (green) 208 green%p3m = .FALSE. 209 green%special_dimension = 0 210 green%radius = 0.0_dp 211 green%slab_size = 0.0_dp 212 green%alpha = 0.0_dp 213 green%method = PERIODIC3D 214 last_greens_fn_id_nr = last_greens_fn_id_nr + 1 215 green%id_nr = last_greens_fn_id_nr 216 green%ref_count = 1 217 green%MT_alpha = 1.0_dp 218 green%MT_rel_cutoff = 1.0_dp 219 green%p3m = .FALSE. 220 green%p3m_order = 0 221 green%p3m_alpha = 0.0_dp 222 green%sr_screening = .FALSE. 223 green%sr_alpha = 1.0_dp 224 green%sr_rc = 0.0_dp 225 226 NULLIFY (green%influence_fn, green%p3m_charge) 227 NULLIFY (green%p3m_coeff, green%p3m_bm2) 228 NULLIFY (green%dct_influence_fn) 229 NULLIFY (green%screen_fn) 230 231 !CPASSERT(cell%orthorhombic) 232 DO i = 1, 3 233 abc(i) = cell_hmat(i, i) 234 END DO 235 dim = COUNT(poisson_params%periodic == 1) 236 237 SELECT CASE (poisson_params%solver) 238 CASE (pw_poisson_periodic) 239 green%method = PERIODIC3D 240 IF (dim /= 3) THEN 241 CPABORT("Illegal combination of periodicity and Poisson solver periodic3d") 242 END IF 243 CASE (pw_poisson_multipole) 244 green%method = MULTIPOLE0D 245 IF (dim /= 0) THEN 246 CPABORT("Illegal combination of periodicity and Poisson solver mulipole0d") 247 END IF 248 CASE (pw_poisson_analytic) 249 SELECT CASE (dim) 250 CASE (0) 251 green%method = ANALYTIC0D 252 green%radius = 0.5_dp*MINVAL(abc) 253 CASE (1) 254 green%method = ANALYTIC1D 255 green%special_dimension = MAXLOC(poisson_params%periodic(1:3), 1) 256 green%radius = MAXVAL(abc(1:3)) 257 DO i = 1, 3 258 IF (i == green%special_dimension) CYCLE 259 green%radius = MIN(green%radius, 0.5_dp*abc(i)) 260 END DO 261 CASE (2) 262 green%method = ANALYTIC2D 263 i = MINLOC(poisson_params%periodic, 1) 264 green%special_dimension = i 265 green%slab_size = abc(i) 266 CASE (3) 267 green%method = PERIODIC3D 268 CASE DEFAULT 269 CPABORT("") 270 END SELECT 271 CASE (pw_poisson_mt) 272 green%MT_rel_cutoff = poisson_params%mt_rel_cutoff 273 green%MT_alpha = poisson_params%mt_alpha 274 green%MT_alpha = green%MT_alpha/MINVAL(abc) 275 SELECT CASE (dim) 276 CASE (0) 277 green%method = MT0D 278 green%radius = 0.5_dp*MINVAL(abc) 279 CASE (1) 280 green%method = MT1D 281 green%special_dimension = MAXLOC(poisson_params%periodic(1:3), 1) 282 green%radius = MAXVAL(abc(1:3)) 283 DO i = 1, 3 284 IF (i == green%special_dimension) CYCLE 285 green%radius = MIN(green%radius, 0.5_dp*abc(i)) 286 END DO 287 CASE (2) 288 green%method = MT2D 289 i = MINLOC(poisson_params%periodic, 1) 290 green%special_dimension = i 291 green%slab_size = abc(i) 292 CASE (3) 293 CPABORT("Illegal combination of periodicity and Poisson solver (MT)") 294 CASE DEFAULT 295 CPABORT("") 296 END SELECT 297 CASE (pw_poisson_implicit) 298 green%method = PS_IMPLICIT 299 CASE DEFAULT 300 CPABORT("An unknown Poisson solver was specified") 301 END SELECT 302 303 ! allocate influence function,... 304 SELECT CASE (green%method) 305 CASE (PERIODIC3D, ANALYTIC2D, ANALYTIC1D, ANALYTIC0D, MT2D, MT1D, MT0D, MULTIPOLE0D, PS_IMPLICIT) 306 CALL pw_pool_create_pw(pw_pool, green%influence_fn, & 307 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 308 309 IF (poisson_params%ewald_type == do_ewald_spme) THEN 310 green%p3m = .TRUE. 311 green%p3m_order = poisson_params%ewald_o_spline 312 green%p3m_alpha = poisson_params%ewald_alpha 313 n = green%p3m_order 314 ALLOCATE (green%p3m_coeff(-(n - 1):n - 1, 0:n - 1)) 315 CALL spme_coeff_calculate(n, green%p3m_coeff) 316 CALL pw_pool_create_pw(pw_pool, green%p3m_charge, use_data=REALDATA1D, & 317 in_space=RECIPROCALSPACE) 318 CALL influence_factor(green) 319 CALL calc_p3m_charge(green) 320 ELSE 321 green%p3m = .FALSE. 322 END IF 323 ! 324 SELECT CASE (green%method) 325 CASE (MT0D, MT1D, MT2D) 326 CALL MTin_create_screen_fn(green%screen_fn, pw_pool=pw_pool, method=green%method, & 327 alpha=green%MT_alpha, & 328 special_dimension=green%special_dimension, slab_size=green%slab_size, & 329 super_ref_pw_grid=mt_super_ref_pw_grid) 330 CASE (PS_IMPLICIT) 331 IF ((poisson_params%ps_implicit_params%boundary_condition .EQ. MIXED_BC) .OR. & 332 (poisson_params%ps_implicit_params%boundary_condition .EQ. NEUMANN_BC)) THEN 333 CALL pw_pool_create(pw_pool_xpndd, pw_grid=dct_pw_grid) 334 CALL pw_pool_create_pw(pw_pool_xpndd, green%dct_influence_fn, & 335 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 336 CALL pw_pool_release(pw_pool_xpndd) 337 END IF 338 END SELECT 339 340 CASE DEFAULT 341 CPABORT("") 342 END SELECT 343 344 ! initialize influence function 345 gf => green%influence_fn 346 grid => green%influence_fn%pw_grid 347 SELECT CASE (green%method) 348 CASE (PERIODIC3D, MULTIPOLE0D) 349 350 DO ig = grid%first_gne0, grid%ngpts_cut_local 351 g2 = grid%gsq(ig) 352 gf%cc(ig) = fourpi/g2 353 END DO 354 IF (grid%have_g0) gf%cc(1) = 0.0_dp 355 356 CASE (ANALYTIC2D) 357 358 iz = green%special_dimension ! iz is the direction with NO PBC 359 zlength = green%slab_size ! zlength is the thickness of the cell 360 DO ig = grid%first_gne0, grid%ngpts_cut_local 361 nz = grid%g_hat(iz, ig) 362 g2 = grid%gsq(ig) 363 g3d = fourpi/g2 364 gg = 0.5_dp*SQRT(g2) 365 gf%cc(ig) = g3d*(1.0_dp - (-1.0_dp)**nz*EXP(-gg*zlength)) 366 END DO 367 IF (grid%have_g0) gf%cc(1) = 0.0_dp 368 369 CASE (ANALYTIC1D) 370 ! see 'ab initio molecular dynamics' table 3.1 371 ! iz is the direction of the PBC ( can be 1,2,3 -> x,y,z ) 372 iz = green%special_dimension 373 ! rlength is the radius of the tube 374 rlength = green%radius 375 DO ig = grid%first_gne0, grid%ngpts_cut_local 376 g2 = grid%gsq(ig) 377 g3d = fourpi/g2 378 gxy = SQRT(MAX(0.0_dp, g2 - grid%g(iz, ig)*grid%g(iz, ig))) 379 gz = ABS(grid%g(iz, ig)) 380 j0g = BESSEL_J0(rlength*gxy) 381 j1g = BESSEL_J1(rlength*gxy) 382 IF (gz > 0) THEN 383 k0g = bessk0(rlength*gz) 384 k1g = bessk1(rlength*gz) 385 ELSE 386 k0g = 0 387 k1g = 0 388 ENDIF 389 gf%cc(ig) = g3d*(1.0_dp + rlength* & 390 (gxy*j1g*k0g - gz*j0g*k1g)) 391 END DO 392 IF (grid%have_g0) gf%cc(1) = 0.0_dp 393 394 CASE (ANALYTIC0D) 395 396 rlength = green%radius ! rlength is the radius of the sphere 397 DO ig = grid%first_gne0, grid%ngpts_cut_local 398 g2 = grid%gsq(ig) 399 gg = SQRT(g2) 400 g3d = fourpi/g2 401 gf%cc(ig) = g3d*(1.0_dp - COS(rlength*gg)) 402 END DO 403 IF (grid%have_g0) & 404 gf%cc(1) = 0.5_dp*fourpi*rlength*rlength 405 406 CASE (MT2D, MT1D, MT0D) 407 408 DO ig = grid%first_gne0, grid%ngpts_cut_local 409 g2 = grid%gsq(ig) 410 g3d = fourpi/g2 411 gf%cc(ig) = g3d+green%screen_fn%cc(ig) 412 END DO 413 IF (grid%have_g0) & 414 gf%cc(1) = green%screen_fn%cc(1) 415 416 CASE (PS_IMPLICIT) 417 418 DO ig = grid%first_gne0, grid%ngpts_cut_local 419 g2 = grid%gsq(ig) 420 gf%cc(ig) = fourpi/g2 421 END DO 422 IF (grid%have_g0) gf%cc(1) = 0.0_dp 423 424 IF (ASSOCIATED(green%dct_influence_fn)) THEN 425 dct_gf => green%dct_influence_fn 426 dct_grid => green%dct_influence_fn%pw_grid 427 428 DO ig = dct_grid%first_gne0, dct_grid%ngpts_cut_local 429 g2 = dct_grid%gsq(ig) 430 dct_gf%cc(ig) = fourpi/g2 431 END DO 432 IF (dct_grid%have_g0) dct_gf%cc(1) = 0.0_dp 433 END IF 434 435 CASE DEFAULT 436 CPABORT("") 437 END SELECT 438 439 END SUBROUTINE pw_green_create 440 441! ************************************************************************************************** 442!> \brief retains the type 443!> \param gftype ... 444!> \author Teodoro Laino 445! ************************************************************************************************** 446 SUBROUTINE pw_green_retain(gftype) 447 TYPE(greens_fn_type), POINTER :: gftype 448 449 CHARACTER(len=*), PARAMETER :: routineN = 'pw_green_retain', & 450 routineP = moduleN//':'//routineN 451 452 CPASSERT(ASSOCIATED(gftype)) 453 CPASSERT(gftype%ref_count > 0) 454 gftype%ref_count = gftype%ref_count + 1 455 END SUBROUTINE pw_green_retain 456 457! ************************************************************************************************** 458!> \brief destroys the type (deallocates data) 459!> \param gftype ... 460!> \param pw_pool ... 461!> \par History 462!> none 463!> \author Joost VandeVondele 464!> Teodoro Laino 465! ************************************************************************************************** 466 SUBROUTINE pw_green_release(gftype, pw_pool) 467 TYPE(greens_fn_type), POINTER :: gftype 468 TYPE(pw_pool_type), OPTIONAL, POINTER :: pw_pool 469 470 CHARACTER(len=*), PARAMETER :: routineN = 'pw_green_release', & 471 routineP = moduleN//':'//routineN 472 473 LOGICAL :: can_give_back 474 475 IF (ASSOCIATED(gftype)) THEN 476 CPASSERT(gftype%ref_count > 0) 477 gftype%ref_count = gftype%ref_count - 1 478 IF (gftype%ref_count == 0) THEN 479 can_give_back = PRESENT(pw_pool) 480 IF (can_give_back) can_give_back = ASSOCIATED(pw_pool) 481 IF (can_give_back) THEN 482 CALL pw_pool_give_back_pw(pw_pool, gftype%influence_fn, & 483 accept_non_compatible=.TRUE.) 484 CALL pw_pool_give_back_pw(pw_pool, gftype%dct_influence_fn, & 485 accept_non_compatible=.TRUE.) 486 CALL pw_pool_give_back_pw(pw_pool, gftype%screen_fn, & 487 accept_non_compatible=.TRUE.) 488 CALL pw_pool_give_back_pw(pw_pool, gftype%p3m_charge, & 489 accept_non_compatible=.TRUE.) 490 ELSE 491 CALL pw_release(gftype%influence_fn) 492 CALL pw_release(gftype%dct_influence_fn) 493 CALL pw_release(gftype%screen_fn) 494 CALL pw_release(gftype%p3m_charge) 495 END IF 496 IF (ASSOCIATED(gftype%p3m_bm2)) & 497 DEALLOCATE (gftype%p3m_bm2) 498 IF (ASSOCIATED(gftype%p3m_coeff)) & 499 DEALLOCATE (gftype%p3m_coeff) 500 DEALLOCATE (gftype) 501 END IF 502 END IF 503 NULLIFY (gftype) 504 END SUBROUTINE pw_green_release 505 506! ************************************************************************************************** 507!> \brief Calculates the influence_factor for the 508!> SPME Green's function in reciprocal space''' 509!> \param gftype ... 510!> \par History 511!> none 512!> \author DH (29-Mar-2001) 513! ************************************************************************************************** 514 SUBROUTINE influence_factor(gftype) 515 TYPE(greens_fn_type), POINTER :: gftype 516 517 CHARACTER(len=*), PARAMETER :: routineN = 'influence_factor', & 518 routineP = moduleN//':'//routineN 519 520 COMPLEX(KIND=dp) :: b_m, exp_m, sum_m 521 INTEGER :: dim, j, k, l, n, pt 522 INTEGER, DIMENSION(3) :: npts 523 INTEGER, DIMENSION(:), POINTER :: lb, ub 524 REAL(KIND=dp) :: l_arg, prod_arg, val 525 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: m_assign 526 527 CPASSERT(ASSOCIATED(gftype)) 528 CPASSERT(gftype%ref_count > 0) 529 n = gftype%p3m_order 530 531 ! calculate the assignment function values 532 533 lb => gftype%influence_fn%pw_grid%bounds(1, :) 534 ub => gftype%influence_fn%pw_grid%bounds(2, :) 535 IF (ASSOCIATED(gftype%p3m_bm2)) THEN 536 IF (LBOUND(gftype%p3m_bm2, 2) /= MINVAL(lb(:)) .OR. & 537 UBOUND(gftype%p3m_bm2, 2) /= MAXVAL(ub(:))) THEN 538 DEALLOCATE (gftype%p3m_bm2) 539 END IF 540 END IF 541 IF (.NOT. ASSOCIATED(gftype%p3m_bm2)) THEN 542 ALLOCATE (gftype%p3m_bm2(3, MINVAL(lb(:)):MAXVAL(ub(:)))) 543 END IF 544 545 ALLOCATE (m_assign(0:n - 2)) 546 m_assign = 0.0_dp 547 DO k = 0, n - 2 548 j = -(n - 1) + 2*k 549 DO l = 0, n - 1 550 l_arg = 0.5_dp**l 551 prod_arg = gftype%p3m_coeff(j, l)*l_arg 552 m_assign(k) = m_assign(k) + prod_arg 553 END DO 554 END DO 555 556 ! calculate the absolute b values 557 558 npts(:) = ub(:) - lb(:) + 1 559 DO dim = 1, 3 560 DO pt = lb(dim), ub(dim) 561 val = twopi*(REAL(pt, KIND=dp)/REAL(npts(dim), KIND=dp)) 562 exp_m = CMPLX(COS(val), -SIN(val), KIND=dp) 563 sum_m = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 564 DO k = 0, n - 2 565 sum_m = sum_m + m_assign(k)*exp_m**k 566 END DO 567 b_m = exp_m**(n - 1)/sum_m 568 gftype%p3m_bm2(dim, pt) = SQRT(REAL(b_m*CONJG(b_m), KIND=dp)) 569 END DO 570 END DO 571 572 DEALLOCATE (m_assign) 573 END SUBROUTINE influence_factor 574 575! ************************************************************************************************** 576!> \brief ... 577!> \param gf ... 578! ************************************************************************************************** 579 SUBROUTINE calc_p3m_charge(gf) 580 581 TYPE(greens_fn_type), POINTER :: gf 582 583 INTEGER :: ig, l, m, n 584 REAL(KIND=dp) :: arg, novol 585 REAL(KIND=dp), DIMENSION(:, :), POINTER :: bm2 586 TYPE(pw_grid_type), POINTER :: grid 587 TYPE(pw_type), POINTER :: pc 588 589 grid => gf%influence_fn%pw_grid 590 591 ! check if charge function is consistent with current box volume 592 593 pc => gf%p3m_charge 594 bm2 => gf%p3m_bm2 595 arg = 1.0_dp/(8.0_dp*gf%p3m_alpha**2) 596 novol = REAL(grid%ngpts, KIND=dp)/grid%vol 597 DO ig = 1, grid%ngpts_cut_local 598 l = grid%g_hat(1, ig) 599 m = grid%g_hat(2, ig) 600 n = grid%g_hat(3, ig) 601 pc%cr(ig) = novol*EXP(-arg*grid%gsq(ig))* & 602 bm2(1, l)*bm2(2, m)*bm2(3, n) 603 END DO 604 605 END SUBROUTINE calc_p3m_charge 606 607! ************************************************************************************************** 608!> \brief Initialize the poisson solver 609!> You should call this just before calling the work routine 610!> pw_poisson_solver 611!> Call pw_poisson_release when you have finished 612!> \param poisson_env ... 613!> \par History 614!> none 615!> \author JGH (12-Mar-2001) 616! ************************************************************************************************** 617 SUBROUTINE pw_poisson_create(poisson_env) 618 619 TYPE(pw_poisson_type), POINTER :: poisson_env 620 621 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_create', & 622 routineP = moduleN//':'//routineN 623 624 CPASSERT(.NOT. ASSOCIATED(poisson_env)) 625 ALLOCATE (poisson_env) 626 last_poisson_id = last_poisson_id + 1 627 poisson_env%id_nr = last_poisson_id 628 poisson_env%ref_count = 1 629 poisson_env%method = pw_poisson_none 630 poisson_env%rebuild = .TRUE. 631 NULLIFY (poisson_env%green_fft, & 632 poisson_env%pw_pools, & 633 poisson_env%mt_super_ref_pw_grid, & 634 poisson_env%wavelet, & 635 poisson_env%implicit_env, & 636 poisson_env%dct_pw_grid, & 637 poisson_env%diel_rs_grid) 638 poisson_env%pw_level = -1 639 poisson_env%ref_count = 1 640 641 END SUBROUTINE pw_poisson_create 642 643! ************************************************************************************************** 644!> \brief retains the pw_poisson_env 645!> \param poisson_env ... 646!> \author fawzi 647! ************************************************************************************************** 648 SUBROUTINE pw_poisson_retain(poisson_env) 649 TYPE(pw_poisson_type), POINTER :: poisson_env 650 651 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_retain', & 652 routineP = moduleN//':'//routineN 653 654 CPASSERT(ASSOCIATED(poisson_env)) 655 CPASSERT(poisson_env%ref_count > 0) 656 poisson_env%ref_count = poisson_env%ref_count + 1 657 END SUBROUTINE pw_poisson_retain 658 659! ************************************************************************************************** 660!> \brief releases the poisson solver 661!> \param poisson_env ... 662!> \par History 663!> none 664!> \author fawzi (11.2002) 665! ************************************************************************************************** 666 SUBROUTINE pw_poisson_release(poisson_env) 667 668 TYPE(pw_poisson_type), POINTER :: poisson_env 669 670 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_release', & 671 routineP = moduleN//':'//routineN 672 673 IF (ASSOCIATED(poisson_env)) THEN 674 CPASSERT(poisson_env%ref_count > 0) 675 poisson_env%ref_count = poisson_env%ref_count - 1 676 IF (poisson_env%ref_count == 0) THEN 677 IF (ASSOCIATED(poisson_env%pw_pools)) THEN 678 CALL pw_pools_dealloc(poisson_env%pw_pools) 679 END IF 680 681 CALL pw_green_release(poisson_env%green_fft) 682 CALL pw_grid_release(poisson_env%mt_super_ref_pw_grid) 683 CALL ps_wavelet_release(poisson_env%wavelet) 684 CALL ps_implicit_release(poisson_env%implicit_env, & 685 poisson_env%parameters%ps_implicit_params) 686 CALL pw_grid_release(poisson_env%dct_pw_grid) 687 CALL rs_grid_release(poisson_env%diel_rs_grid) 688 DEALLOCATE (poisson_env) 689 690 END IF 691 END IF 692 NULLIFY (poisson_env) 693 694 END SUBROUTINE pw_poisson_release 695 696! ************************************************************************************************** 697!> \brief Calculates the coefficients for the charge assignment function 698!> \param n ... 699!> \param coeff ... 700!> \par History 701!> none 702!> \author DG (29-Mar-2001) 703! ************************************************************************************************** 704 SUBROUTINE spme_coeff_calculate(n, coeff) 705 706 INTEGER, INTENT(IN) :: n 707 REAL(KIND=dp), DIMENSION(-(n-1):n-1, 0:n-1), & 708 INTENT(OUT) :: coeff 709 710 INTEGER :: i, j, l, m 711 REAL(KIND=dp) :: b 712 REAL(KIND=dp), DIMENSION(n, -n:n, 0:n-1) :: a 713 714 a = 0.0_dp 715 a(1, 0, 0) = 1.0_dp 716 717 DO i = 2, n 718 m = i - 1 719 DO j = -m, m, 2 720 DO l = 0, m - 1 721 b = (a(m, j - 1, l) + & 722 REAL((-1)**l, KIND=dp)*a(m, j + 1, l))/ & 723 REAL((l + 1)*2**(l + 1), KIND=dp) 724 a(i, j, 0) = a(i, j, 0) + b 725 END DO 726 DO l = 0, m - 1 727 a(i, j, l + 1) = (a(m, j + 1, l) - & 728 a(m, j - 1, l))/REAL(l + 1, KIND=dp) 729 END DO 730 END DO 731 END DO 732 733 coeff = 0.0_dp 734 DO i = 0, n - 1 735 DO j = -(n - 1), n - 1, 2 736 coeff(j, i) = a(n, j, i) 737 END DO 738 END DO 739 740 END SUBROUTINE spme_coeff_calculate 741 742END MODULE pw_poisson_types 743