1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculate the electrostatic energy by the Smooth Particle Ewald method 8!> \par History 9!> JGH (03-May-2001) : first correctly working version 10!> \author JGH (21-Mar-2001) 11! ************************************************************************************************** 12MODULE spme 13 14 USE atomic_kind_types, ONLY: get_atomic_kind 15 USE atprop_types, ONLY: atprop_type 16 USE bibliography, ONLY: Essmann1995,& 17 cite_reference 18 USE cell_types, ONLY: cell_type,& 19 real_to_scaled 20 USE dgs, ONLY: dg_sum_patch,& 21 dg_sum_patch_force_1d,& 22 dg_sum_patch_force_3d 23 USE ewald_environment_types, ONLY: ewald_env_get,& 24 ewald_environment_type 25 USE ewald_pw_types, ONLY: ewald_pw_get,& 26 ewald_pw_type 27 USE kinds, ONLY: dp 28 USE mathconstants, ONLY: fourpi 29 USE particle_types, ONLY: particle_type 30 USE pme_tools, ONLY: get_center,& 31 set_list 32 USE pw_grid_types, ONLY: pw_grid_type 33 USE pw_grids, ONLY: get_pw_grid_info 34 USE pw_methods, ONLY: pw_copy,& 35 pw_derive,& 36 pw_integral_a2b,& 37 pw_transfer 38 USE pw_poisson_methods, ONLY: pw_poisson_rebuild,& 39 pw_poisson_solve 40 USE pw_poisson_types, ONLY: greens_fn_type,& 41 pw_poisson_type 42 USE pw_pool_types, ONLY: pw_pool_create_pw,& 43 pw_pool_give_back_pw,& 44 pw_pool_type 45 USE pw_types, ONLY: COMPLEXDATA1D,& 46 REALDATA3D,& 47 REALSPACE,& 48 RECIPROCALSPACE,& 49 pw_p_type,& 50 pw_type 51 USE realspace_grid_types, ONLY: & 52 pw2rs, realspace_grid_desc_type, realspace_grid_p_type, realspace_grid_type, rs2pw, & 53 rs_grid_create, rs_grid_release, rs_grid_set_box, rs_grid_zero, rs_pw_transfer 54 USE shell_potential_types, ONLY: shell_kind_type 55#include "./base/base_uses.f90" 56 57 IMPLICIT NONE 58 59 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'spme' 60 61 PRIVATE 62 PUBLIC :: spme_evaluate, spme_potential, spme_forces, get_patch 63 64 INTERFACE get_patch 65 MODULE PROCEDURE get_patch_a, get_patch_b 66 END INTERFACE 67 68CONTAINS 69 70! ************************************************************************************************** 71!> \brief ... 72!> \param ewald_env ... 73!> \param ewald_pw ... 74!> \param box ... 75!> \param particle_set ... 76!> \param fg_coulomb ... 77!> \param vg_coulomb ... 78!> \param pv_g ... 79!> \param shell_particle_set ... 80!> \param core_particle_set ... 81!> \param fgshell_coulomb ... 82!> \param fgcore_coulomb ... 83!> \param use_virial ... 84!> \param charges ... 85!> \param atprop ... 86!> \par History 87!> JGH (03-May-2001) : SPME with charge definition 88!> \author JGH (21-Mar-2001) 89! ************************************************************************************************** 90 SUBROUTINE spme_evaluate(ewald_env, ewald_pw, box, particle_set, & 91 fg_coulomb, vg_coulomb, pv_g, shell_particle_set, core_particle_set, & 92 fgshell_coulomb, fgcore_coulomb, use_virial, charges, atprop) 93 94 TYPE(ewald_environment_type), POINTER :: ewald_env 95 TYPE(ewald_pw_type), POINTER :: ewald_pw 96 TYPE(cell_type), POINTER :: box 97 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set 98 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: fg_coulomb 99 REAL(KIND=dp), INTENT(OUT) :: vg_coulomb 100 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: pv_g 101 TYPE(particle_type), DIMENSION(:), OPTIONAL, & 102 POINTER :: shell_particle_set, core_particle_set 103 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), & 104 OPTIONAL :: fgshell_coulomb, fgcore_coulomb 105 LOGICAL, INTENT(IN) :: use_virial 106 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges 107 TYPE(atprop_type), POINTER :: atprop 108 109 CHARACTER(len=*), PARAMETER :: routineN = 'spme_evaluate', routineP = moduleN//':'//routineN 110 111 INTEGER :: group, handle, i, ig, ipart, j, n, & 112 ncore, npart, nshell, o_spline, p1, & 113 p1_shell 114 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: center, core_center, shell_center 115 INTEGER, DIMENSION(3) :: nd, npts 116 LOGICAL :: do_shell 117 REAL(KIND=dp) :: alpha, dvols, fat1, ffa, ffb 118 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rhos 119 REAL(KIND=dp), DIMENSION(3) :: fat 120 REAL(KIND=dp), DIMENSION(3, 3) :: f_stress, h_stress 121 TYPE(greens_fn_type), POINTER :: green 122 TYPE(pw_grid_type), POINTER :: grid_spme 123 TYPE(pw_p_type), DIMENSION(3) :: dphi_g 124 TYPE(pw_poisson_type), POINTER :: poisson_env 125 TYPE(pw_pool_type), POINTER :: pw_pool 126 TYPE(pw_type), POINTER :: phi_g, rhob_g, rhob_r 127 TYPE(realspace_grid_desc_type), POINTER :: rs_desc 128 TYPE(realspace_grid_p_type), DIMENSION(:), POINTER :: drpot 129 TYPE(realspace_grid_type), POINTER :: rden, rpot 130 131 NULLIFY (drpot, grid_spme, green, poisson_env, phi_g, rhob_g, rhob_r, pw_pool, rden, rpot) 132 CALL timeset(routineN, handle) 133 CALL cite_reference(Essmann1995) 134 135 !-------------- INITIALISATION --------------------- 136 CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, & 137 group=group) 138 CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, & 139 poisson_env=poisson_env) 140 CALL pw_poisson_rebuild(poisson_env) 141 green => poisson_env%green_fft 142 grid_spme => pw_pool%pw_grid 143 144 npart = SIZE(particle_set) 145 146 CALL get_pw_grid_info(grid_spme, npts=npts, dvol=dvols) 147 148 n = o_spline 149 ALLOCATE (rhos(n, n, n)) 150 CALL rs_grid_create(rden, rs_desc) 151 CALL rs_grid_set_box(grid_spme, rs=rden) 152 CALL rs_grid_zero(rden) 153 154 ALLOCATE (center(3, npart)) 155 CALL get_center(particle_set, box, center, npts, n) 156 IF (PRESENT(shell_particle_set) .AND. (PRESENT(core_particle_set))) THEN 157 CPASSERT(ASSOCIATED(shell_particle_set)) 158 CPASSERT(ASSOCIATED(core_particle_set)) 159 nshell = SIZE(shell_particle_set) 160 ncore = SIZE(core_particle_set) 161 CPASSERT(nshell == ncore) 162 ALLOCATE (shell_center(3, nshell)) 163 CALL get_center(shell_particle_set, box, shell_center, npts, n) 164 ALLOCATE (core_center(3, nshell)) 165 CALL get_center(core_particle_set, box, core_center, npts, n) 166 END IF 167 168 !-------------- DENSITY CALCULATION ---------------- 169 ipart = 0 170 ! Particles 171 DO 172 CALL set_list(particle_set, npart, center, p1, rden, ipart) 173 IF (p1 == 0) EXIT 174 175 do_shell = (particle_set(p1)%shell_index /= 0) 176 IF (do_shell) CYCLE 177 ! calculate function on small boxes 178 CALL get_patch(particle_set, box, green, npts, p1, rhos, is_core=.FALSE., & 179 is_shell=.FALSE., unit_charge=.FALSE., charges=charges) 180 181 ! add boxes to real space grid (big box) 182 CALL dg_sum_patch(rden, rhos, center(:, p1)) 183 END DO 184 ! Shell-Model 185 IF (PRESENT(shell_particle_set) .AND. PRESENT(core_particle_set)) THEN 186 ipart = 0 187 DO 188 ! calculate function on small boxes 189 CALL set_list(shell_particle_set, nshell, shell_center, p1_shell, & 190 rden, ipart) 191 IF (p1_shell == 0) EXIT 192 CALL get_patch(shell_particle_set, box, green, npts, p1_shell, rhos, & 193 is_core=.FALSE., is_shell=.TRUE., unit_charge=.FALSE.) 194 195 ! add boxes to real space grid (big box) 196 CALL dg_sum_patch(rden, rhos, shell_center(:, p1_shell)) 197 END DO 198 ipart = 0 199 DO 200 ! calculate function on small boxes 201 CALL set_list(core_particle_set, nshell, core_center, p1_shell, & 202 rden, ipart) 203 IF (p1_shell == 0) EXIT 204 CALL get_patch(core_particle_set, box, green, npts, p1_shell, rhos, & 205 is_core=.TRUE., is_shell=.FALSE., unit_charge=.FALSE.) 206 207 ! add boxes to real space grid (big box) 208 CALL dg_sum_patch(rden, rhos, core_center(:, p1_shell)) 209 END DO 210 END IF 211 !----------- END OF DENSITY CALCULATION ------------- 212 213 CALL pw_pool_create_pw(pw_pool, rhob_r, use_data=REALDATA3D, & 214 in_space=REALSPACE) 215 CALL rs_pw_transfer(rden, rhob_r, rs2pw) 216 ! transform density to G space and add charge function 217 CALL pw_pool_create_pw(pw_pool, rhob_g, & 218 use_data=COMPLEXDATA1D, & 219 in_space=RECIPROCALSPACE) 220 CALL pw_transfer(rhob_r, rhob_g) 221 ! update charge function 222 rhob_g%cc = rhob_g%cc*green%p3m_charge%cr 223 224 !-------------- ELECTROSTATIC CALCULATION ----------- 225 ! allocate intermediate arrays 226 DO i = 1, 3 227 NULLIFY (dphi_g(i)%pw) 228 CALL pw_pool_create_pw(pw_pool, dphi_g(i)%pw, & 229 use_data=COMPLEXDATA1D, & 230 in_space=RECIPROCALSPACE) 231 END DO 232 CALL pw_pool_create_pw(pw_pool, phi_g, & 233 use_data=COMPLEXDATA1D, & 234 in_space=RECIPROCALSPACE) 235 CALL pw_poisson_solve(poisson_env, rhob_g, vg_coulomb, phi_g, dphi_g, & 236 h_stress) 237 !---------- END OF ELECTROSTATIC CALCULATION -------- 238 239 ! Atomic Energy and Stress 240 IF (atprop%energy .OR. atprop%stress) THEN 241 CALL rs_grid_create(rpot, rs_desc) 242 CALL rs_grid_set_box(grid_spme, rs=rpot) 243 phi_g%cc = phi_g%cc*green%p3m_charge%cr 244 CALL pw_transfer(phi_g, rhob_r) 245 CALL rs_pw_transfer(rpot, rhob_r, pw2rs) 246 ipart = 0 247 DO 248 CALL set_list(particle_set, npart, center, p1, rden, ipart) 249 IF (p1 == 0) EXIT 250 do_shell = (particle_set(p1)%shell_index /= 0) 251 IF (do_shell) CYCLE 252 ! calculate function on small boxes 253 CALL get_patch(particle_set, box, green, grid_spme%npts, p1, rhos, is_core=.FALSE., & 254 is_shell=.FALSE., unit_charge=.FALSE., charges=charges) 255 ! integrate box and potential 256 CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fat1) 257 IF (atprop%energy) THEN 258 atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols 259 END IF 260 IF (atprop%stress) THEN 261 atprop%atstress(1, 1, p1) = atprop%atstress(1, 1, p1) + 0.5_dp*fat1*dvols 262 atprop%atstress(2, 2, p1) = atprop%atstress(2, 2, p1) + 0.5_dp*fat1*dvols 263 atprop%atstress(3, 3, p1) = atprop%atstress(3, 3, p1) + 0.5_dp*fat1*dvols 264 END IF 265 END DO 266 ! Core-shell model 267 IF (PRESENT(shell_particle_set) .AND. PRESENT(core_particle_set)) THEN 268 ipart = 0 269 DO 270 CALL set_list(shell_particle_set, nshell, shell_center, p1_shell, rden, ipart) 271 IF (p1_shell == 0) EXIT 272 CALL get_patch(shell_particle_set, box, green, npts, p1_shell, rhos, & 273 is_core=.FALSE., is_shell=.TRUE., unit_charge=.FALSE.) 274 CALL dg_sum_patch_force_1d(rpot, rhos, shell_center(:, p1_shell), fat1) 275 p1 = shell_particle_set(p1_shell)%atom_index 276 IF (atprop%energy) THEN 277 atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols 278 END IF 279 END DO 280 ipart = 0 281 DO 282 CALL set_list(core_particle_set, nshell, core_center, p1_shell, rden, ipart) 283 IF (p1_shell == 0) EXIT 284 CALL get_patch(core_particle_set, box, green, npts, p1_shell, rhos, & 285 is_core=.TRUE., is_shell=.FALSE., unit_charge=.FALSE.) 286 CALL dg_sum_patch_force_1d(rpot, rhos, core_center(:, p1_shell), fat1) 287 p1 = core_particle_set(p1_shell)%atom_index 288 IF (atprop%energy) THEN 289 atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols 290 END IF 291 END DO 292 END IF 293 IF (atprop%stress) THEN 294 ffa = (0.5_dp/alpha)**2 295 ffb = 1.0_dp/fourpi 296 DO i = 1, 3 297 DO ig = grid_spme%first_gne0, grid_spme%ngpts_cut_local 298 phi_g%cc(ig) = ffb*dphi_g(i)%pw%cc(ig)*(ffa*grid_spme%gsq(ig) + 1.0_dp) 299 phi_g%cc(ig) = phi_g%cc(ig)*poisson_env%green_fft%influence_fn%cc(ig) 300 END DO 301 IF (grid_spme%have_g0) phi_g%cc(1) = 0.0_dp 302 DO j = 1, i 303 nd = 0 304 nd(j) = 1 305 CALL pw_copy(phi_g, rhob_g) 306 CALL pw_derive(rhob_g, nd) 307 rhob_g%cc = rhob_g%cc*green%p3m_charge%cr 308 CALL pw_transfer(rhob_g, rhob_r) 309 CALL rs_pw_transfer(rpot, rhob_r, pw2rs) 310 311 ipart = 0 312 DO 313 CALL set_list(particle_set, npart, center, p1, rden, ipart) 314 IF (p1 == 0) EXIT 315 ! calculate function on small boxes 316 CALL get_patch(particle_set, box, green, grid_spme%npts, p1, rhos, & 317 is_core=.FALSE., is_shell=.FALSE., unit_charge=.FALSE., charges=charges) 318 ! integrate box and potential 319 CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fat1) 320 atprop%atstress(i, j, p1) = atprop%atstress(i, j, p1) + fat1*dvols 321 IF (i /= j) atprop%atstress(j, i, p1) = atprop%atstress(j, i, p1) + fat1*dvols 322 END DO 323 324 END DO 325 END DO 326 END IF 327 CALL rs_grid_release(rpot) 328 END IF 329 330 CALL pw_pool_give_back_pw(pw_pool, phi_g) 331 CALL pw_pool_give_back_pw(pw_pool, rhob_g) 332 333 !------------- STRESS TENSOR CALCULATION ------------ 334 IF (use_virial) THEN 335 DO i = 1, 3 336 DO j = i, 3 337 f_stress(i, j) = pw_integral_a2b(dphi_g(i)%pw, dphi_g(j)%pw) 338 f_stress(j, i) = f_stress(i, j) 339 END DO 340 END DO 341 ffa = (1.0_dp/fourpi)*(0.5_dp/alpha)**2 342 f_stress = -ffa*f_stress 343 pv_g = h_stress + f_stress 344 END IF 345 !--------END OF STRESS TENSOR CALCULATION ----------- 346 ! move derivative of potential to real space grid and 347 ! multiply by charge function in g-space 348 ALLOCATE (drpot(1:3)) 349 DO i = 1, 3 350 CALL rs_grid_create(drpot(i)%rs_grid, rs_desc) 351 CALL rs_grid_set_box(grid_spme, rs=drpot(i)%rs_grid) 352 dphi_g(i)%pw%cc = dphi_g(i)%pw%cc*green%p3m_charge%cr 353 CALL pw_transfer(dphi_g(i)%pw, rhob_r) 354 CALL pw_pool_give_back_pw(pw_pool, dphi_g(i)%pw) 355 CALL rs_pw_transfer(drpot(i)%rs_grid, rhob_r, pw2rs) 356 END DO 357 358 CALL pw_pool_give_back_pw(pw_pool, rhob_r) 359 !----------------- FORCE CALCULATION ---------------- 360 ! initialize the forces 361 fg_coulomb = 0.0_dp 362 ! Particles 363 ipart = 0 364 DO 365 CALL set_list(particle_set, npart, center, p1, rden, ipart) 366 IF (p1 == 0) EXIT 367 368 do_shell = (particle_set(p1)%shell_index /= 0) 369 IF (do_shell) CYCLE 370 ! calculate function on small boxes 371 CALL get_patch(particle_set, box, green, grid_spme%npts, p1, rhos, is_core=.FALSE., & 372 is_shell=.FALSE., unit_charge=.FALSE., charges=charges) 373 374 ! add boxes to real space grid (big box) 375 CALL dg_sum_patch_force_3d(drpot, rhos, center(:, p1), fat) 376 fg_coulomb(1, p1) = fg_coulomb(1, p1) - fat(1)*dvols 377 fg_coulomb(2, p1) = fg_coulomb(2, p1) - fat(2)*dvols 378 fg_coulomb(3, p1) = fg_coulomb(3, p1) - fat(3)*dvols 379 END DO 380 ! Shell-Model 381 IF (PRESENT(shell_particle_set) .AND. (PRESENT(core_particle_set))) THEN 382 IF (PRESENT(fgshell_coulomb)) THEN 383 ipart = 0 384 fgshell_coulomb = 0.0_dp 385 DO 386 ! calculate function on small boxes 387 CALL set_list(shell_particle_set, nshell, shell_center, p1_shell, & 388 rden, ipart) 389 IF (p1_shell == 0) EXIT 390 391 CALL get_patch(shell_particle_set, box, green, grid_spme%npts, & 392 p1_shell, rhos, is_core=.FALSE., is_shell=.TRUE., unit_charge=.FALSE.) 393 394 ! add boxes to real space grid (big box) 395 CALL dg_sum_patch_force_3d(drpot, rhos, shell_center(:, p1_shell), fat) 396 fgshell_coulomb(1, p1_shell) = fgshell_coulomb(1, p1_shell) - fat(1)*dvols 397 fgshell_coulomb(2, p1_shell) = fgshell_coulomb(2, p1_shell) - fat(2)*dvols 398 fgshell_coulomb(3, p1_shell) = fgshell_coulomb(3, p1_shell) - fat(3)*dvols 399 400 END DO 401 END IF 402 IF (PRESENT(fgcore_coulomb)) THEN 403 ipart = 0 404 fgcore_coulomb = 0.0_dp 405 DO 406 ! calculate function on small boxes 407 CALL set_list(core_particle_set, nshell, core_center, p1_shell, & 408 rden, ipart) 409 IF (p1_shell == 0) EXIT 410 411 CALL get_patch(core_particle_set, box, green, grid_spme%npts, & 412 p1_shell, rhos, is_core=.TRUE., is_shell=.FALSE., unit_charge=.FALSE.) 413 414 ! add boxes to real space grid (big box) 415 CALL dg_sum_patch_force_3d(drpot, rhos, core_center(:, p1_shell), fat) 416 fgcore_coulomb(1, p1_shell) = fgcore_coulomb(1, p1_shell) - fat(1)*dvols 417 fgcore_coulomb(2, p1_shell) = fgcore_coulomb(2, p1_shell) - fat(2)*dvols 418 fgcore_coulomb(3, p1_shell) = fgcore_coulomb(3, p1_shell) - fat(3)*dvols 419 END DO 420 END IF 421 422 END IF 423 !--------------END OF FORCE CALCULATION ------------- 424 425 !------------------CLEANING UP ---------------------- 426 CALL rs_grid_release(rden) 427 IF (ASSOCIATED(drpot)) THEN 428 DO i = 1, 3 429 CALL rs_grid_release(drpot(i)%rs_grid) 430 END DO 431 DEALLOCATE (drpot) 432 END IF 433 434 DEALLOCATE (rhos) 435 DEALLOCATE (center) 436 IF (ALLOCATED(shell_center)) THEN 437 DEALLOCATE (shell_center) 438 END IF 439 IF (ALLOCATED(core_center)) THEN 440 DEALLOCATE (core_center) 441 END IF 442 CALL timestop(handle) 443 444 END SUBROUTINE spme_evaluate 445 446! ************************************************************************************************** 447!> \brief Calculate the electrostatic potential from particles A (charge A) at 448!> positions of particles B 449!> \param ewald_env ... 450!> \param ewald_pw ... 451!> \param box ... 452!> \param particle_set_a ... 453!> \param charges_a ... 454!> \param particle_set_b ... 455!> \param potential ... 456!> \par History 457!> JGH (23-Oct-2014) : adapted from SPME evaluate 458!> \author JGH (23-Oct-2014) 459! ************************************************************************************************** 460 SUBROUTINE spme_potential(ewald_env, ewald_pw, box, particle_set_a, charges_a, & 461 particle_set_b, potential) 462 463 TYPE(ewald_environment_type), POINTER :: ewald_env 464 TYPE(ewald_pw_type), POINTER :: ewald_pw 465 TYPE(cell_type), POINTER :: box 466 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set_a 467 REAL(KIND=dp), DIMENSION(:), POINTER :: charges_a 468 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set_b 469 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: potential 470 471 CHARACTER(len=*), PARAMETER :: routineN = 'spme_potential', routineP = moduleN//':'//routineN 472 473 INTEGER :: group, handle, ipart, n, npart_a, & 474 npart_b, o_spline, p1 475 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: center 476 INTEGER, DIMENSION(3) :: npts 477 REAL(KIND=dp) :: alpha, dvols, fat1 478 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rhos 479 TYPE(greens_fn_type), POINTER :: green 480 TYPE(pw_grid_type), POINTER :: grid_spme 481 TYPE(pw_poisson_type), POINTER :: poisson_env 482 TYPE(pw_pool_type), POINTER :: pw_pool 483 TYPE(pw_type), POINTER :: phi_g, rhob_g, rhob_r 484 TYPE(realspace_grid_desc_type), POINTER :: rs_desc 485 TYPE(realspace_grid_type), POINTER :: rden, rpot 486 487 NULLIFY (grid_spme, green, poisson_env, phi_g, rhob_g, rhob_r, pw_pool, & 488 rden, rpot) 489 CALL timeset(routineN, handle) 490 CALL cite_reference(Essmann1995) 491 492 !-------------- INITIALISATION --------------------- 493 CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group) 494 CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, poisson_env=poisson_env) 495 CALL pw_poisson_rebuild(poisson_env) 496 green => poisson_env%green_fft 497 grid_spme => pw_pool%pw_grid 498 499 npart_a = SIZE(particle_set_a) 500 npart_b = SIZE(particle_set_b) 501 502 CALL get_pw_grid_info(grid_spme, npts=npts, dvol=dvols) 503 504 n = o_spline 505 ALLOCATE (rhos(n, n, n)) 506 CALL rs_grid_create(rden, rs_desc) 507 CALL rs_grid_set_box(grid_spme, rs=rden) 508 CALL rs_grid_zero(rden) 509 510 ALLOCATE (center(3, npart_a)) 511 CALL get_center(particle_set_a, box, center, npts, n) 512 513 !-------------- DENSITY CALCULATION ---------------- 514 ipart = 0 515 ! Particles 516 DO 517 CALL set_list(particle_set_a, npart_a, center, p1, rden, ipart) 518 IF (p1 == 0) EXIT 519 520 ! calculate function on small boxes 521 CALL get_patch(particle_set_a, box, green, npts, p1, rhos, charges_a) 522 523 ! add boxes to real space grid (big box) 524 CALL dg_sum_patch(rden, rhos, center(:, p1)) 525 END DO 526 DEALLOCATE (center) 527 !----------- END OF DENSITY CALCULATION ------------- 528 529 CALL pw_pool_create_pw(pw_pool, rhob_r, use_data=REALDATA3D, & 530 in_space=REALSPACE) 531 CALL rs_pw_transfer(rden, rhob_r, rs2pw) 532 ! transform density to G space and add charge function 533 CALL pw_pool_create_pw(pw_pool, rhob_g, & 534 use_data=COMPLEXDATA1D, & 535 in_space=RECIPROCALSPACE) 536 CALL pw_transfer(rhob_r, rhob_g) 537 ! update charge function 538 rhob_g%cc = rhob_g%cc*green%p3m_charge%cr 539 540 !-------------- ELECTROSTATIC CALCULATION ----------- 541 ! allocate intermediate arrays 542 CALL pw_pool_create_pw(pw_pool, phi_g, & 543 use_data=COMPLEXDATA1D, & 544 in_space=RECIPROCALSPACE) 545 CALL pw_poisson_solve(poisson_env, density=rhob_g, vhartree=phi_g) 546 !---------- END OF ELECTROSTATIC CALCULATION -------- 547 548 !-------------- POTENTAIL AT ATOMIC POSITIONS ------- 549 ALLOCATE (center(3, npart_b)) 550 CALL get_center(particle_set_b, box, center, npts, n) 551 rpot => rden 552 phi_g%cc = phi_g%cc*green%p3m_charge%cr 553 CALL pw_transfer(phi_g, rhob_r) 554 CALL rs_pw_transfer(rpot, rhob_r, pw2rs) 555 ipart = 0 556 DO 557 CALL set_list(particle_set_b, npart_b, center, p1, rden, ipart) 558 IF (p1 == 0) EXIT 559 ! calculate function on small boxes 560 CALL get_patch(particle_set_b, box, green, grid_spme%npts, p1, rhos, & 561 is_core=.FALSE., is_shell=.FALSE., unit_charge=.TRUE.) 562 ! integrate box and potential 563 CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fat1) 564 potential(p1) = potential(p1) + fat1*dvols 565 END DO 566 567 !------------------CLEANING UP ---------------------- 568 CALL pw_pool_give_back_pw(pw_pool, phi_g) 569 CALL pw_pool_give_back_pw(pw_pool, rhob_g) 570 CALL pw_pool_give_back_pw(pw_pool, rhob_r) 571 CALL rs_grid_release(rden) 572 573 DEALLOCATE (rhos) 574 DEALLOCATE (center) 575 CALL timestop(handle) 576 577 END SUBROUTINE spme_potential 578 579! ************************************************************************************************** 580!> \brief Calculate the forces on particles B for the electrostatic interaction 581!> betrween particles A and B 582!> \param ewald_env ... 583!> \param ewald_pw ... 584!> \param box ... 585!> \param particle_set_a ... 586!> \param charges_a ... 587!> \param particle_set_b ... 588!> \param charges_b ... 589!> \param forces_b ... 590!> \par History 591!> JGH (27-Oct-2014) : adapted from SPME evaluate 592!> \author JGH (23-Oct-2014) 593! ************************************************************************************************** 594 SUBROUTINE spme_forces(ewald_env, ewald_pw, box, particle_set_a, charges_a, & 595 particle_set_b, charges_b, forces_b) 596 597 TYPE(ewald_environment_type), POINTER :: ewald_env 598 TYPE(ewald_pw_type), POINTER :: ewald_pw 599 TYPE(cell_type), POINTER :: box 600 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set_a 601 REAL(KIND=dp), DIMENSION(:), POINTER :: charges_a 602 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set_b 603 REAL(KIND=dp), DIMENSION(:), POINTER :: charges_b 604 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: forces_b 605 606 CHARACTER(len=*), PARAMETER :: routineN = 'spme_forces', routineP = moduleN//':'//routineN 607 608 INTEGER :: group, handle, i, ipart, n, npart_a, & 609 npart_b, o_spline, p1 610 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: center 611 INTEGER, DIMENSION(3) :: npts 612 REAL(KIND=dp) :: alpha, dvols 613 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rhos 614 REAL(KIND=dp), DIMENSION(3) :: fat 615 TYPE(greens_fn_type), POINTER :: green 616 TYPE(pw_grid_type), POINTER :: grid_spme 617 TYPE(pw_p_type), DIMENSION(3) :: dphi_g 618 TYPE(pw_poisson_type), POINTER :: poisson_env 619 TYPE(pw_pool_type), POINTER :: pw_pool 620 TYPE(pw_type), POINTER :: phi_g, rhob_g, rhob_r 621 TYPE(realspace_grid_desc_type), POINTER :: rs_desc 622 TYPE(realspace_grid_p_type), DIMENSION(:), POINTER :: drpot 623 TYPE(realspace_grid_type), POINTER :: rden 624 625 NULLIFY (drpot, grid_spme, green, poisson_env, phi_g, rhob_g, rhob_r, & 626 pw_pool, rden) 627 CALL timeset(routineN, handle) 628 CALL cite_reference(Essmann1995) 629 630 !-------------- INITIALISATION --------------------- 631 CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group) 632 CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, poisson_env=poisson_env) 633 CALL pw_poisson_rebuild(poisson_env) 634 green => poisson_env%green_fft 635 grid_spme => pw_pool%pw_grid 636 637 npart_a = SIZE(particle_set_a) 638 npart_b = SIZE(particle_set_b) 639 640 CALL get_pw_grid_info(grid_spme, npts=npts, dvol=dvols) 641 642 n = o_spline 643 ALLOCATE (rhos(n, n, n)) 644 CALL rs_grid_create(rden, rs_desc) 645 CALL rs_grid_set_box(grid_spme, rs=rden) 646 CALL rs_grid_zero(rden) 647 648 ALLOCATE (center(3, npart_a)) 649 CALL get_center(particle_set_a, box, center, npts, n) 650 651 !-------------- DENSITY CALCULATION ---------------- 652 ipart = 0 653 ! Particles 654 DO 655 CALL set_list(particle_set_a, npart_a, center, p1, rden, ipart) 656 IF (p1 == 0) EXIT 657 658 ! calculate function on small boxes 659 CALL get_patch(particle_set_a, box, green, npts, p1, rhos, charges_a) 660 661 ! add boxes to real space grid (big box) 662 CALL dg_sum_patch(rden, rhos, center(:, p1)) 663 END DO 664 DEALLOCATE (center) 665 !----------- END OF DENSITY CALCULATION ------------- 666 667 CALL pw_pool_create_pw(pw_pool, rhob_r, use_data=REALDATA3D, & 668 in_space=REALSPACE) 669 CALL rs_pw_transfer(rden, rhob_r, rs2pw) 670 ! transform density to G space and add charge function 671 CALL pw_pool_create_pw(pw_pool, rhob_g, & 672 use_data=COMPLEXDATA1D, & 673 in_space=RECIPROCALSPACE) 674 CALL pw_transfer(rhob_r, rhob_g) 675 ! update charge function 676 rhob_g%cc = rhob_g%cc*green%p3m_charge%cr 677 678 !-------------- ELECTROSTATIC CALCULATION ----------- 679 ! allocate intermediate arrays 680 DO i = 1, 3 681 NULLIFY (dphi_g(i)%pw) 682 CALL pw_pool_create_pw(pw_pool, dphi_g(i)%pw, & 683 use_data=COMPLEXDATA1D, & 684 in_space=RECIPROCALSPACE) 685 END DO 686 CALL pw_pool_create_pw(pw_pool, phi_g, & 687 use_data=COMPLEXDATA1D, & 688 in_space=RECIPROCALSPACE) 689 CALL pw_poisson_solve(poisson_env, density=rhob_g, vhartree=phi_g, & 690 dvhartree=dphi_g) 691 !---------- END OF ELECTROSTATIC CALCULATION -------- 692 ! move derivative of potential to real space grid and 693 ! multiply by charge function in g-space 694 ALLOCATE (drpot(1:3)) 695 DO i = 1, 3 696 CALL rs_grid_create(drpot(i)%rs_grid, rs_desc) 697 CALL rs_grid_set_box(grid_spme, rs=drpot(i)%rs_grid) 698 dphi_g(i)%pw%cc = dphi_g(i)%pw%cc*green%p3m_charge%cr 699 CALL pw_transfer(dphi_g(i)%pw, rhob_r) 700 CALL pw_pool_give_back_pw(pw_pool, dphi_g(i)%pw) 701 CALL rs_pw_transfer(drpot(i)%rs_grid, rhob_r, pw2rs) 702 END DO 703 !----------------- FORCE CALCULATION ---------------- 704 ALLOCATE (center(3, npart_b)) 705 CALL get_center(particle_set_b, box, center, npts, n) 706 ipart = 0 707 DO 708 CALL set_list(particle_set_b, npart_b, center, p1, rden, ipart) 709 IF (p1 == 0) EXIT 710 ! calculate function on small boxes 711 CALL get_patch(particle_set_b, box, green, grid_spme%npts, p1, rhos, & 712 is_core=.FALSE., is_shell=.FALSE., unit_charge=.FALSE., charges=charges_b) 713 ! add boxes to real space grid (big box) 714 CALL dg_sum_patch_force_3d(drpot, rhos, center(:, p1), fat) 715 forces_b(1, p1) = forces_b(1, p1) - fat(1)*dvols 716 forces_b(2, p1) = forces_b(2, p1) - fat(2)*dvols 717 forces_b(3, p1) = forces_b(3, p1) - fat(3)*dvols 718 END DO 719 !------------------CLEANING UP ---------------------- 720 IF (ASSOCIATED(drpot)) THEN 721 DO i = 1, 3 722 CALL rs_grid_release(drpot(i)%rs_grid) 723 END DO 724 DEALLOCATE (drpot) 725 END IF 726 CALL pw_pool_give_back_pw(pw_pool, phi_g) 727 CALL pw_pool_give_back_pw(pw_pool, rhob_g) 728 CALL pw_pool_give_back_pw(pw_pool, rhob_r) 729 CALL rs_grid_release(rden) 730 731 DEALLOCATE (rhos) 732 DEALLOCATE (center) 733 CALL timestop(handle) 734 735 END SUBROUTINE spme_forces 736 737! ************************************************************************************************** 738!> \brief Calculates local density in a small box 739!> \param part ... 740!> \param box ... 741!> \param green ... 742!> \param npts ... 743!> \param p ... 744!> \param rhos ... 745!> \param is_core ... 746!> \param is_shell ... 747!> \param unit_charge ... 748!> \param charges ... 749!> \par History 750!> none 751!> \author JGH (21-Mar-2001) 752! ************************************************************************************************** 753 SUBROUTINE get_patch_a(part, box, green, npts, p, rhos, is_core, is_shell, & 754 unit_charge, charges) 755 756 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: part 757 TYPE(cell_type), POINTER :: box 758 TYPE(greens_fn_type), POINTER :: green 759 INTEGER, DIMENSION(3), INTENT(IN) :: npts 760 INTEGER, INTENT(IN) :: p 761 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: rhos 762 LOGICAL, INTENT(IN) :: is_core, is_shell, unit_charge 763 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges 764 765 CHARACTER(len=*), PARAMETER :: routineN = 'get_patch_a', routineP = moduleN//':'//routineN 766 767 INTEGER :: nbox 768 LOGICAL :: use_charge_array 769 REAL(KIND=dp) :: q 770 REAL(KIND=dp), DIMENSION(3) :: delta, r 771 TYPE(shell_kind_type), POINTER :: shell 772 773 NULLIFY (shell) 774 use_charge_array = PRESENT(charges) 775 IF (use_charge_array) use_charge_array = ASSOCIATED(charges) 776 IF (is_core .AND. is_shell) THEN 777 CPABORT("Shell-model: cannot be core and shell simultaneously") 778 END IF 779 780 nbox = SIZE(rhos, 1) 781 r = part(p)%r 782 q = 1.0_dp 783 IF (.NOT. unit_charge) THEN 784 IF (is_core) THEN 785 CALL get_atomic_kind(atomic_kind=part(p)%atomic_kind, shell=shell) 786 q = shell%charge_core 787 ELSE IF (is_shell) THEN 788 CALL get_atomic_kind(atomic_kind=part(p)%atomic_kind, shell=shell) 789 q = shell%charge_shell 790 ELSE 791 CALL get_atomic_kind(atomic_kind=part(p)%atomic_kind, qeff=q) 792 END IF 793 IF (use_charge_array) q = charges(p) 794 END IF 795 CALL get_delta(box, r, npts, delta, nbox) 796 CALL spme_get_patch(rhos, nbox, delta, q, green%p3m_coeff) 797 798 END SUBROUTINE get_patch_a 799 800! ************************************************************************************************** 801!> \brief Calculates local density in a small box 802!> \param part ... 803!> \param box ... 804!> \param green ... 805!> \param npts ... 806!> \param p ... 807!> \param rhos ... 808!> \param charges ... 809!> \par History 810!> none 811!> \author JGH (21-Mar-2001) 812! ************************************************************************************************** 813 SUBROUTINE get_patch_b(part, box, green, npts, p, rhos, charges) 814 815 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: part 816 TYPE(cell_type), POINTER :: box 817 TYPE(greens_fn_type), POINTER :: green 818 INTEGER, DIMENSION(3), INTENT(IN) :: npts 819 INTEGER, INTENT(IN) :: p 820 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: rhos 821 REAL(KIND=dp), DIMENSION(:), POINTER :: charges 822 823 CHARACTER(len=*), PARAMETER :: routineN = 'get_patch_b', routineP = moduleN//':'//routineN 824 825 INTEGER :: nbox 826 REAL(KIND=dp) :: q 827 REAL(KIND=dp), DIMENSION(3) :: delta, r 828 829 nbox = SIZE(rhos, 1) 830 r = part(p)%r 831 q = charges(p) 832 CALL get_delta(box, r, npts, delta, nbox) 833 CALL spme_get_patch(rhos, nbox, delta, q, green%p3m_coeff) 834 835 END SUBROUTINE get_patch_b 836 837! ************************************************************************************************** 838!> \brief Calculates SPME charge assignment 839!> \param rhos ... 840!> \param n ... 841!> \param delta ... 842!> \param q ... 843!> \param coeff ... 844!> \par History 845!> DG (29-Mar-2001) : code implemented 846!> \author JGH (22-Mar-2001) 847! ************************************************************************************************** 848 SUBROUTINE spme_get_patch(rhos, n, delta, q, coeff) 849 850 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: rhos 851 INTEGER, INTENT(IN) :: n 852 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: delta 853 REAL(KIND=dp), INTENT(IN) :: q 854 REAL(KIND=dp), DIMENSION(-(n-1):n-1, 0:n-1), & 855 INTENT(IN) :: coeff 856 857 CHARACTER(len=*), PARAMETER :: routineN = 'spme_get_patch', routineP = moduleN//':'//routineN 858 INTEGER, PARAMETER :: nmax = 12 859 860 INTEGER :: i, i1, i2, i3, j, l 861 REAL(KIND=dp) :: r2, r3 862 REAL(KIND=dp), DIMENSION(3, -nmax:nmax) :: w_assign 863 REAL(KIND=dp), DIMENSION(3, 0:nmax-1) :: deltal 864 REAL(KIND=dp), DIMENSION(3, 1:nmax) :: f_assign 865 866 IF (n > nmax) THEN 867 CPABORT("nmax value too small") 868 END IF 869 ! calculate the assignment function values and 870 ! the charges on the grid (small box) 871 872 deltal(1, 0) = 1.0_dp 873 deltal(2, 0) = 1.0_dp 874 deltal(3, 0) = 1.0_dp 875 DO l = 1, n - 1 876 deltal(1, l) = deltal(1, l - 1)*delta(1) 877 deltal(2, l) = deltal(2, l - 1)*delta(2) 878 deltal(3, l) = deltal(3, l - 1)*delta(3) 879 END DO 880 881 w_assign = 0.0_dp 882 DO j = -(n - 1), n - 1, 2 883 DO l = 0, n - 1 884 w_assign(1, j) = w_assign(1, j) + coeff(j, l)*deltal(1, l) 885 w_assign(2, j) = w_assign(2, j) + coeff(j, l)*deltal(2, l) 886 w_assign(3, j) = w_assign(3, j) + coeff(j, l)*deltal(3, l) 887 END DO 888 END DO 889 DO i = 1, n 890 j = n + 1 - 2*i 891 f_assign(1, i) = w_assign(1, j) 892 f_assign(2, i) = w_assign(2, j) 893 f_assign(3, i) = w_assign(3, j) 894 END DO 895 896 DO i3 = 1, n 897 r3 = q*f_assign(3, i3) 898 DO i2 = 1, n 899 r2 = r3*f_assign(2, i2) 900 DO i1 = 1, n 901 rhos(i1, i2, i3) = r2*f_assign(1, i1) 902 END DO 903 END DO 904 END DO 905 906 END SUBROUTINE spme_get_patch 907 908! ************************************************************************************************** 909!> \brief ... 910!> \param box ... 911!> \param r ... 912!> \param npts ... 913!> \param delta ... 914!> \param n ... 915! ************************************************************************************************** 916 SUBROUTINE get_delta(box, r, npts, delta, n) 917 918 TYPE(cell_type), POINTER :: box 919 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r 920 INTEGER, DIMENSION(3), INTENT(IN) :: npts 921 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: delta 922 INTEGER, INTENT(IN) :: n 923 924 CHARACTER(len=*), PARAMETER :: routineN = 'get_delta', routineP = moduleN//':'//routineN 925 926 INTEGER :: mp 927 INTEGER, DIMENSION(3) :: center 928 REAL(KIND=dp) :: rmp 929 REAL(KIND=dp), DIMENSION(3) :: ca, grid_i, s 930 931 mp = MAXVAL(npts(:)) 932 rmp = REAL(mp, KIND=dp) 933 ! compute the scaled coordinate of atomi 934 CALL real_to_scaled(s, r, box) 935 s = s - REAL(NINT(s), KIND=dp) 936 937 ! find the continuous ``grid'' point 938 grid_i(1:3) = REAL(npts(1:3), KIND=dp)*s(1:3) 939 940 ! find the closest grid point 941 942 IF (MOD(n, 2) == 0) THEN 943 center(:) = INT(grid_i(:) + rmp) - mp 944 ca(:) = REAL(center(:), KIND=dp) + 0.5_dp 945 ELSE 946 center(:) = NINT(grid_i(:)) 947 ca(:) = REAL(center(:), KIND=dp) 948 END IF 949 950 ! find the distance vector 951 delta(:) = grid_i(:) - ca(:) 952 953 END SUBROUTINE get_delta 954 955END MODULE spme 956 957