1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> JGH (21-Mar-2001) : Complete rewrite 9!> \author CJM and APSI 10! ************************************************************************************************** 11MODULE pme 12 13 USE atomic_kind_types, ONLY: atomic_kind_type,& 14 get_atomic_kind 15 USE atprop_types, ONLY: atprop_type 16 USE bibliography, ONLY: cite_reference,& 17 darden1993 18 USE cell_types, ONLY: cell_type 19 USE dg_rho0_types, ONLY: dg_rho0_type 20 USE dg_types, ONLY: dg_get,& 21 dg_type 22 USE dgs, ONLY: dg_get_patch,& 23 dg_get_strucfac,& 24 dg_sum_patch,& 25 dg_sum_patch_force_1d,& 26 dg_sum_patch_force_3d 27 USE ewald_environment_types, ONLY: ewald_env_get,& 28 ewald_environment_type 29 USE ewald_pw_types, ONLY: ewald_pw_get,& 30 ewald_pw_type 31 USE kinds, ONLY: dp 32 USE mathconstants, ONLY: fourpi 33 USE particle_types, ONLY: particle_type 34 USE pme_tools, ONLY: get_center,& 35 set_list 36 USE pw_grid_types, ONLY: pw_grid_type 37 USE pw_methods, ONLY: pw_copy,& 38 pw_derive,& 39 pw_integral_a2b,& 40 pw_transfer 41 USE pw_poisson_methods, ONLY: pw_poisson_solve 42 USE pw_poisson_types, ONLY: pw_poisson_type 43 USE pw_pool_types, ONLY: pw_pool_create_pw,& 44 pw_pool_give_back_pw,& 45 pw_pool_type 46 USE pw_types, ONLY: COMPLEXDATA1D,& 47 REALDATA3D,& 48 REALSPACE,& 49 RECIPROCALSPACE,& 50 pw_p_type,& 51 pw_type 52 USE realspace_grid_types, ONLY: & 53 pw2rs, realspace_grid_desc_type, realspace_grid_p_type, realspace_grid_type, rs2pw, & 54 rs_grid_create, rs_grid_release, rs_grid_set_box, rs_grid_zero, rs_pw_transfer 55 USE shell_potential_types, ONLY: shell_kind_type 56 USE structure_factor_types, ONLY: structure_factor_type 57 USE structure_factors, ONLY: structure_factor_allocate,& 58 structure_factor_deallocate,& 59 structure_factor_init 60#include "./base/base_uses.f90" 61 62 IMPLICIT NONE 63 64 PRIVATE 65 PUBLIC :: pme_evaluate 66 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pme' 67 68CONTAINS 69 70! ************************************************************************************************** 71!> \brief ... 72!> \param ewald_env ... 73!> \param ewald_pw ... 74!> \param box ... 75!> \param particle_set ... 76!> \param vg_coulomb ... 77!> \param fg_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 (15-Mar-2001) : New electrostatic calculation and pressure tensor 88!> JGH (21-Mar-2001) : Complete rewrite 89!> JGH (21-Mar-2001) : Introduced real space density type for future 90!> parallelisation 91!> \author CJM and APSI 92! ************************************************************************************************** 93 SUBROUTINE pme_evaluate(ewald_env, ewald_pw, box, particle_set, vg_coulomb, & 94 fg_coulomb, pv_g, shell_particle_set, core_particle_set, & 95 fgshell_coulomb, fgcore_coulomb, use_virial, charges, atprop) 96 TYPE(ewald_environment_type), POINTER :: ewald_env 97 TYPE(ewald_pw_type), POINTER :: ewald_pw 98 TYPE(cell_type), POINTER :: box 99 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 100 REAL(KIND=dp), INTENT(OUT) :: vg_coulomb 101 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), & 102 OPTIONAL :: fg_coulomb, pv_g 103 TYPE(particle_type), DIMENSION(:), OPTIONAL, & 104 POINTER :: shell_particle_set, core_particle_set 105 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), & 106 OPTIONAL :: fgshell_coulomb, fgcore_coulomb 107 LOGICAL, INTENT(IN) :: use_virial 108 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges 109 TYPE(atprop_type), POINTER :: atprop 110 111 CHARACTER(LEN=*), PARAMETER :: routineN = 'pme_evaluate', routineP = moduleN//':'//routineN 112 113 INTEGER :: group, handle, i, ig, ipart, j, nd(3), & 114 npart, nshell, p1, p2 115 LOGICAL :: is1_core, is2_core 116 REAL(KIND=dp) :: alpha, dvols, fat1, ffa, ffb 117 REAL(KIND=dp), DIMENSION(3) :: fat 118 REAL(KIND=dp), DIMENSION(3, 3) :: f_stress, h_stress 119 TYPE(dg_rho0_type), POINTER :: dg_rho0 120 TYPE(dg_type), POINTER :: dg 121 TYPE(pw_grid_type), POINTER :: grid_b, grid_s 122 TYPE(pw_p_type) :: rhos1, rhos2 123 TYPE(pw_p_type), DIMENSION(3) :: dphi_g 124 TYPE(pw_poisson_type), POINTER :: poisson_env 125 TYPE(pw_pool_type), POINTER :: pw_big_pool, pw_small_pool 126 TYPE(pw_type), POINTER :: phi_g, phi_r, 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 TYPE(structure_factor_type) :: exp_igr 131 132 CALL timeset(routineN, handle) 133 NULLIFY (poisson_env, rden, drpot) 134 CALL cite_reference(Darden1993) 135 CALL ewald_env_get(ewald_env, alpha=alpha, group=group) 136 CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_big_pool, & 137 pw_small_pool=pw_small_pool, rs_desc=rs_desc, & 138 poisson_env=poisson_env, dg=dg) 139 140 grid_b => pw_big_pool%pw_grid 141 grid_s => pw_small_pool%pw_grid 142 143 CALL dg_get(dg, dg_rho0=dg_rho0) 144 145 npart = SIZE(particle_set) 146 147 CALL structure_factor_init(exp_igr) 148 149 IF (PRESENT(shell_particle_set)) THEN 150 CPASSERT(ASSOCIATED(shell_particle_set)) 151 CPASSERT(ASSOCIATED(core_particle_set)) 152 nshell = SIZE(shell_particle_set) 153 CALL structure_factor_allocate(grid_s%bounds, npart, exp_igr, & 154 allocate_centre=.TRUE., allocate_shell_e=.TRUE., & 155 allocate_shell_centre=.TRUE., nshell=nshell) 156 157 ELSE 158 CALL structure_factor_allocate(grid_s%bounds, npart, exp_igr, & 159 allocate_centre=.TRUE.) 160 END IF 161 162 CALL pw_pool_create_pw(pw_small_pool, rhos1%pw, & 163 use_data=REALDATA3D) 164 CALL pw_pool_create_pw(pw_small_pool, rhos2%pw, & 165 use_data=REALDATA3D) 166 167 CALL rs_grid_create(rden, rs_desc) 168 CALL rs_grid_set_box(grid_b, rs=rden) 169 CALL rs_grid_zero(rden) 170 171 IF (rden%desc%parallel .AND. rden%desc%distributed) THEN 172 CALL get_center(particle_set, box, exp_igr%centre, grid_b%npts, 1) 173 END IF 174 IF (PRESENT(shell_particle_set) .AND. rden%desc%parallel .AND. rden%desc%distributed) THEN 175 CALL get_center(shell_particle_set, box, exp_igr%shell_centre, grid_b%npts, 1) 176 CALL get_center(core_particle_set, box, exp_igr%core_centre, grid_b%npts, 1) 177 END IF 178 179 !-------------- DENSITY CALCULATION ---------------- 180 181 ipart = 0 182 DO 183 184 CALL set_list(particle_set, npart, exp_igr%centre, p1, rden, ipart, exp_igr%core_centre) 185 CALL set_list(particle_set, npart, exp_igr%centre, p2, rden, ipart, exp_igr%core_centre) 186 IF (p1 == 0 .AND. p2 == 0) EXIT 187 188 is1_core = (particle_set(p1)%shell_index /= 0) 189 IF (p2 /= 0) THEN 190 is2_core = (particle_set(p2)%shell_index /= 0) 191 ELSE 192 is2_core = .FALSE. 193 END IF 194 195 ! calculate function on small boxes (we use double packing in FFT) 196 IF (is1_core .OR. is2_core) THEN 197 CALL get_patch(dg, particle_set, exp_igr, box, p1, p2, grid_b, grid_s, & 198 rhos1, rhos2, is1_core=is1_core, is2_core=is2_core, & 199 core_particle_set=core_particle_set, charges=charges) 200 201 ! add boxes to real space grid (big box) 202 IF (is1_core) THEN 203 CALL dg_sum_patch(rden, rhos1, exp_igr%core_centre(:, particle_set(p1)%shell_index)) 204 ELSE 205 CALL dg_sum_patch(rden, rhos1, exp_igr%centre(:, p1)) 206 END IF 207 IF (p2 /= 0 .AND. is2_core) THEN 208 CALL dg_sum_patch(rden, rhos2, exp_igr%core_centre(:, particle_set(p2)%shell_index)) 209 ELSE IF (p2 /= 0) THEN 210 CALL dg_sum_patch(rden, rhos2, exp_igr%centre(:, p2)) 211 END IF 212 ELSE 213 CALL get_patch(dg, particle_set, exp_igr, box, p1, p2, grid_b, grid_s, & 214 rhos1, rhos2, charges=charges) 215 ! add boxes to real space grid (big box) 216 CALL dg_sum_patch(rden, rhos1, exp_igr%centre(:, p1)) 217 IF (p2 /= 0) CALL dg_sum_patch(rden, rhos2, exp_igr%centre(:, p2)) 218 END IF 219 220 END DO 221 IF (PRESENT(shell_particle_set)) THEN 222 ipart = 0 223 DO 224 CALL set_list(shell_particle_set, nshell, exp_igr%shell_centre, p1, rpot, ipart) 225 CALL set_list(shell_particle_set, nshell, exp_igr%shell_centre, p2, rpot, ipart) 226 IF (p1 == 0 .AND. p2 == 0) EXIT 227 ! calculate function on small boxes (we use double packing in FFT) 228 CALL get_patch(dg, shell_particle_set, exp_igr, box, p1, p2, grid_b, grid_s, & 229 rhos1, rhos2, is1_shell=.TRUE., is2_shell=.TRUE., charges=charges) 230 ! add boxes to real space grid (big box) 231 CALL dg_sum_patch(rpot, rhos1, exp_igr%shell_centre(:, p1)) 232 IF (p2 /= 0) CALL dg_sum_patch(rpot, rhos2, exp_igr%shell_centre(:, p2)) 233 END DO 234 END IF 235 236 CALL pw_pool_create_pw(pw_big_pool, rhob_r, use_data=REALDATA3D, in_space=REALSPACE) 237 CALL rs_pw_transfer(rden, rhob_r, rs2pw) 238 239 !-------------- ELECTROSTATIC CALCULATION ----------- 240 241 ! allocate intermediate arrays 242 DO i = 1, 3 243 NULLIFY (dphi_g(i)%pw) 244 CALL pw_pool_create_pw(pw_big_pool, dphi_g(i)%pw, & 245 use_data=COMPLEXDATA1D, & 246 in_space=RECIPROCALSPACE) 247 END DO 248 CALL pw_pool_create_pw(pw_big_pool, phi_r, & 249 use_data=REALDATA3D, & 250 in_space=REALSPACE) 251 252 CALL pw_poisson_solve(poisson_env, rhob_r, vg_coulomb, phi_r, dphi_g, h_stress) 253 254 ! atomic energies 255 IF (atprop%energy .OR. atprop%stress) THEN 256 dvols = rhos1%pw%pw_grid%dvol 257 CALL rs_grid_create(rpot, rs_desc) 258 CALL rs_pw_transfer(rpot, phi_r, pw2rs) 259 ipart = 0 260 DO 261 CALL set_list(particle_set, npart, exp_igr%centre, p1, rden, ipart, exp_igr%core_centre) 262 CALL set_list(particle_set, npart, exp_igr%centre, p2, rden, ipart, exp_igr%core_centre) 263 IF (p1 == 0 .AND. p2 == 0) EXIT 264 ! integrate box and potential 265 CALL get_patch(dg, particle_set, exp_igr, box, p1, p2, grid_b, grid_s, & 266 rhos1, rhos2, charges=charges) 267 ! add boxes to real space grid (big box) 268 CALL dg_sum_patch_force_1d(rpot, rhos1, exp_igr%centre(:, p1), fat1) 269 IF (atprop%energy) THEN 270 atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols 271 END IF 272 IF (atprop%stress) THEN 273 atprop%atstress(1, 1, p1) = atprop%atstress(1, 1, p1) + 0.5_dp*fat1*dvols 274 atprop%atstress(2, 2, p1) = atprop%atstress(2, 2, p1) + 0.5_dp*fat1*dvols 275 atprop%atstress(3, 3, p1) = atprop%atstress(3, 3, p1) + 0.5_dp*fat1*dvols 276 END IF 277 IF (p2 /= 0) THEN 278 CALL dg_sum_patch_force_1d(rpot, rhos2, exp_igr%centre(:, p2), fat1) 279 IF (atprop%energy) THEN 280 atprop%atener(p2) = atprop%atener(p2) + 0.5_dp*fat1*dvols 281 END IF 282 IF (atprop%stress) THEN 283 atprop%atstress(1, 1, p2) = atprop%atstress(1, 1, p2) + 0.5_dp*fat1*dvols 284 atprop%atstress(2, 2, p2) = atprop%atstress(2, 2, p2) + 0.5_dp*fat1*dvols 285 atprop%atstress(3, 3, p2) = atprop%atstress(3, 3, p2) + 0.5_dp*fat1*dvols 286 END IF 287 END IF 288 END DO 289 IF (atprop%stress) THEN 290 CALL pw_pool_create_pw(pw_big_pool, phi_g, & 291 use_data=COMPLEXDATA1D, & 292 in_space=RECIPROCALSPACE) 293 CALL pw_pool_create_pw(pw_big_pool, rhob_g, & 294 use_data=COMPLEXDATA1D, & 295 in_space=RECIPROCALSPACE) 296 ffa = (0.5_dp/dg_rho0%zet(1))**2 297 ffb = 1.0_dp/fourpi 298 DO i = 1, 3 299 DO ig = grid_b%first_gne0, grid_b%ngpts_cut_local 300 phi_g%cc(ig) = ffb*dphi_g(i)%pw%cc(ig)*(ffa*grid_b%gsq(ig) + 1.0_dp) 301 phi_g%cc(ig) = phi_g%cc(ig)*poisson_env%green_fft%influence_fn%cc(ig) 302 END DO 303 IF (grid_b%have_g0) phi_g%cc(1) = 0.0_dp 304 DO j = 1, i 305 CALL pw_copy(phi_g, rhob_g) 306 nd = 0 307 nd(j) = 1 308 CALL pw_derive(rhob_g, nd) 309 CALL pw_transfer(rhob_g, rhob_r) 310 CALL rs_pw_transfer(rpot, rhob_r, pw2rs) 311 312 ipart = 0 313 DO 314 CALL set_list(particle_set, npart, exp_igr%centre, p1, rden, ipart, exp_igr%core_centre) 315 CALL set_list(particle_set, npart, exp_igr%centre, p2, rden, ipart, exp_igr%core_centre) 316 IF (p1 == 0 .AND. p2 == 0) EXIT 317 ! integrate box and potential 318 CALL get_patch(dg, particle_set, exp_igr, box, p1, p2, grid_b, grid_s, & 319 rhos1, rhos2, charges=charges) 320 ! add boxes to real space grid (big box) 321 CALL dg_sum_patch_force_1d(rpot, rhos1, exp_igr%centre(:, p1), fat1) 322 atprop%atstress(i, j, p1) = atprop%atstress(i, j, p1) + fat1*dvols 323 IF (i /= j) atprop%atstress(j, i, p1) = atprop%atstress(j, i, p1) + fat1*dvols 324 IF (p2 /= 0) THEN 325 CALL dg_sum_patch_force_1d(rpot, rhos2, exp_igr%centre(:, p2), fat1) 326 atprop%atstress(i, j, p2) = atprop%atstress(i, j, p2) + fat1*dvols 327 IF (i /= j) atprop%atstress(j, i, p2) = atprop%atstress(j, i, p2) + fat1*dvols 328 END IF 329 END DO 330 END DO 331 END DO 332 CALL pw_pool_give_back_pw(pw_big_pool, phi_g) 333 CALL pw_pool_give_back_pw(pw_big_pool, rhob_g) 334 END IF 335 CALL rs_grid_release(rpot) 336 END IF 337 338 CALL pw_pool_give_back_pw(pw_big_pool, phi_r) 339 340 !---------- END OF ELECTROSTATIC CALCULATION -------- 341 342 !------------- STRESS TENSOR CALCULATION ------------ 343 344 IF ((use_virial) .AND. (PRESENT(pv_g))) THEN 345 DO i = 1, 3 346 DO j = i, 3 347 f_stress(i, j) = pw_integral_a2b(dphi_g(i)%pw, dphi_g(j)%pw) 348 f_stress(j, i) = f_stress(i, j) 349 END DO 350 END DO 351 ffa = (1.0_dp/fourpi)*(0.5_dp/dg_rho0%zet(1))**2 352 f_stress = -ffa*f_stress 353 pv_g = h_stress + f_stress 354 END IF 355 356 !--------END OF STRESS TENSOR CALCULATION ----------- 357 358 ALLOCATE (drpot(1:3)) 359 DO i = 1, 3 360 CALL rs_grid_create(drpot(i)%rs_grid, rs_desc) 361 CALL rs_grid_set_box(grid_b, rs=drpot(i)%rs_grid) 362 CALL pw_transfer(dphi_g(i)%pw, rhob_r) 363 CALL pw_pool_give_back_pw(pw_big_pool, dphi_g(i)%pw) 364 CALL rs_pw_transfer(drpot(i)%rs_grid, rhob_r, pw2rs) 365 END DO 366 367 CALL pw_pool_give_back_pw(pw_big_pool, rhob_r) 368 369 !----------------- FORCE CALCULATION ---------------- 370 371 ! initialize the forces 372 IF (PRESENT(fg_coulomb)) THEN 373 fg_coulomb = 0.0_dp 374 dvols = rhos1%pw%pw_grid%dvol 375 376 ipart = 0 377 DO 378 379 CALL set_list(particle_set, npart, exp_igr%centre, p1, rden, ipart, exp_igr%core_centre) 380 CALL set_list(particle_set, npart, exp_igr%centre, p2, rden, ipart, exp_igr%core_centre) 381 IF (p1 == 0 .AND. p2 == 0) EXIT 382 383 is1_core = (particle_set(p1)%shell_index /= 0) 384 IF (p2 /= 0) THEN 385 is2_core = (particle_set(p2)%shell_index /= 0) 386 ELSE 387 is2_core = .FALSE. 388 END IF 389 390 ! calculate function on small boxes (we use double packing in FFT) 391 392 CALL get_patch_again(dg, particle_set, exp_igr, p1, p2, rhos1, rhos2, & 393 is1_core=is1_core, is2_core=is2_core, charges=charges) 394 395 ! sum boxes on real space grids (big box) 396 IF (is1_core) THEN 397 CALL dg_sum_patch_force_3d(drpot, rhos1, & 398 exp_igr%core_centre(:, particle_set(p1)%shell_index), fat) 399 fgcore_coulomb(1, particle_set(p1)%shell_index) = & 400 fgcore_coulomb(1, particle_set(p1)%shell_index) - fat(1)*dvols 401 fgcore_coulomb(2, particle_set(p1)%shell_index) = & 402 fgcore_coulomb(2, particle_set(p1)%shell_index) - fat(2)*dvols 403 fgcore_coulomb(3, particle_set(p1)%shell_index) = & 404 fgcore_coulomb(3, particle_set(p1)%shell_index) - fat(3)*dvols 405 ELSE 406 CALL dg_sum_patch_force_3d(drpot, rhos1, exp_igr%centre(:, p1), fat) 407 fg_coulomb(1, p1) = fg_coulomb(1, p1) - fat(1)*dvols 408 fg_coulomb(2, p1) = fg_coulomb(2, p1) - fat(2)*dvols 409 fg_coulomb(3, p1) = fg_coulomb(3, p1) - fat(3)*dvols 410 END IF 411 IF (p2 /= 0 .AND. is2_core) THEN 412 CALL dg_sum_patch_force_3d(drpot, rhos1, & 413 exp_igr%core_centre(:, particle_set(p2)%shell_index), fat) 414 fgcore_coulomb(1, particle_set(p2)%shell_index) = & 415 fgcore_coulomb(1, particle_set(p2)%shell_index) - fat(1)*dvols 416 fgcore_coulomb(2, particle_set(p2)%shell_index) = & 417 fgcore_coulomb(2, particle_set(p2)%shell_index) - fat(2)*dvols 418 fgcore_coulomb(3, particle_set(p2)%shell_index) = & 419 fgcore_coulomb(3, particle_set(p2)%shell_index) - fat(3)*dvols 420 ELSEIF (p2 /= 0) THEN 421 CALL dg_sum_patch_force_3d(drpot, rhos2, exp_igr%centre(:, p2), fat) 422 fg_coulomb(1, p2) = fg_coulomb(1, p2) - fat(1)*dvols 423 fg_coulomb(2, p2) = fg_coulomb(2, p2) - fat(2)*dvols 424 fg_coulomb(3, p2) = fg_coulomb(3, p2) - fat(3)*dvols 425 END IF 426 427 END DO 428 ENDIF 429 IF (PRESENT(fgshell_coulomb)) THEN 430 fgshell_coulomb = 0.0_dp 431 dvols = rhos1%pw%pw_grid%dvol 432 433 ipart = 0 434 DO 435 CALL set_list(shell_particle_set, nshell, exp_igr%shell_centre, p1, rden, ipart) 436 CALL set_list(shell_particle_set, nshell, exp_igr%shell_centre, p2, rden, ipart) 437 IF (p1 == 0 .AND. p2 == 0) EXIT 438 439 ! calculate function on small boxes (we use double packing in FFT) 440 CALL get_patch_again(dg, shell_particle_set, exp_igr, p1, p2, rhos1, rhos2, & 441 is1_shell=.TRUE., is2_shell=.TRUE., charges=charges) 442 443 ! sum boxes on real space grids (big box) 444 CALL dg_sum_patch_force_3d(drpot, rhos1, exp_igr%shell_centre(:, p1), fat) 445 fgshell_coulomb(1, p1) = fgshell_coulomb(1, p1) - fat(1)*dvols 446 fgshell_coulomb(2, p1) = fgshell_coulomb(2, p1) - fat(2)*dvols 447 fgshell_coulomb(3, p1) = fgshell_coulomb(3, p1) - fat(3)*dvols 448 IF (p2 /= 0) THEN 449 CALL dg_sum_patch_force_3d(drpot, rhos2, exp_igr%shell_centre(:, p2), fat) 450 fgshell_coulomb(1, p2) = fgshell_coulomb(1, p2) - fat(1)*dvols 451 fgshell_coulomb(2, p2) = fgshell_coulomb(2, p2) - fat(2)*dvols 452 fgshell_coulomb(3, p2) = fgshell_coulomb(3, p2) - fat(3)*dvols 453 END IF 454 END DO 455 456 END IF 457 !--------------END OF FORCE CALCULATION ------------- 458 459 !------------------CLEANING UP ---------------------- 460 461 CALL rs_grid_release(rden) 462 IF (ASSOCIATED(drpot)) THEN 463 DO i = 1, 3 464 CALL rs_grid_release(drpot(i)%rs_grid) 465 END DO 466 DEALLOCATE (drpot) 467 END IF 468 469 CALL pw_pool_give_back_pw(pw_small_pool, rhos1%pw) 470 CALL pw_pool_give_back_pw(pw_small_pool, rhos2%pw) 471 CALL structure_factor_deallocate(exp_igr) 472 473 CALL timestop(handle) 474 475 END SUBROUTINE pme_evaluate 476 477! ************************************************************************************************** 478!> \brief Calculates local density in a small box 479!> \param dg ... 480!> \param particle_set ... 481!> \param exp_igr ... 482!> \param box ... 483!> \param p1 ... 484!> \param p2 ... 485!> \param grid_b ... 486!> \param grid_s ... 487!> \param rhos1 ... 488!> \param rhos2 ... 489!> \param is1_core ... 490!> \param is2_core ... 491!> \param is1_shell ... 492!> \param is2_shell ... 493!> \param core_particle_set ... 494!> \param charges ... 495!> \par History 496!> JGH (23-Mar-2001) : Switch to integer from particle list pointers 497!> \author JGH (21-Mar-2001) 498! ************************************************************************************************** 499 SUBROUTINE get_patch(dg, particle_set, exp_igr, box, p1, p2, & 500 grid_b, grid_s, rhos1, rhos2, is1_core, is2_core, is1_shell, & 501 is2_shell, core_particle_set, charges) 502 503 TYPE(dg_type), POINTER :: dg 504 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 505 TYPE(structure_factor_type) :: exp_igr 506 TYPE(cell_type), INTENT(IN) :: box 507 INTEGER, INTENT(IN) :: p1, p2 508 TYPE(pw_grid_type), INTENT(IN) :: grid_b, grid_s 509 TYPE(pw_p_type) :: rhos1, rhos2 510 LOGICAL, OPTIONAL :: is1_core, is2_core, is1_shell, is2_shell 511 TYPE(particle_type), DIMENSION(:), OPTIONAL, & 512 POINTER :: core_particle_set 513 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges 514 515 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_patch', routineP = moduleN//':'//routineN 516 517 COMPLEX(KIND=dp), DIMENSION(:), POINTER :: ex1, ex2, ey1, ey2, ez1, ez2 518 INTEGER, DIMENSION(:), POINTER :: center1, center2 519 LOGICAL :: my_is1_core, my_is1_shell, my_is2_core, & 520 my_is2_shell, use_charge_array 521 REAL(KIND=dp) :: q1, q2 522 REAL(KIND=dp), DIMENSION(3) :: r1, r2 523 TYPE(atomic_kind_type), POINTER :: atomic_kind 524 TYPE(dg_rho0_type), POINTER :: dg_rho0 525 TYPE(pw_p_type), POINTER :: rho0 526 TYPE(shell_kind_type), POINTER :: shell 527 528 NULLIFY (shell) 529 use_charge_array = PRESENT(charges) 530 IF (use_charge_array) use_charge_array = ASSOCIATED(charges) 531 my_is1_core = .FALSE. 532 my_is2_core = .FALSE. 533 IF (PRESENT(is1_core)) my_is1_core = is1_core 534 IF (PRESENT(is2_core)) my_is2_core = is2_core 535 IF (my_is1_core .OR. my_is2_core) THEN 536 CPASSERT(PRESENT(core_particle_set)) 537 END IF 538 my_is1_shell = .FALSE. 539 my_is2_shell = .FALSE. 540 IF (PRESENT(is1_shell)) my_is1_shell = is1_shell 541 IF (PRESENT(is2_shell)) my_is2_shell = is2_shell 542 IF (my_is1_core .AND. my_is1_shell) THEN 543 CPABORT("Shell-model: cannot be core and shell simultaneously") 544 END IF 545 546 CALL dg_get(dg, dg_rho0=dg_rho0) 547 rho0 => dg_rho0%density 548 549 IF (my_is1_core) THEN 550 r1 = core_particle_set(particle_set(p1)%shell_index)%r 551 ELSE 552 r1 = particle_set(p1)%r 553 END IF 554 atomic_kind => particle_set(p1)%atomic_kind 555 IF (my_is1_core) THEN 556 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell) 557 q1 = shell%charge_core 558 ELSE IF (my_is1_shell) THEN 559 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell) 560 q1 = shell%charge_shell 561 ELSE 562 CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q1) 563 END IF 564 IF (use_charge_array) q1 = charges(p1) 565 566 IF (my_is1_shell) THEN 567 center1 => exp_igr%shell_centre(:, p1) 568 ex1 => exp_igr%shell_ex(:, p1) 569 ey1 => exp_igr%shell_ey(:, p1) 570 ez1 => exp_igr%shell_ez(:, p1) 571 ELSEIF (my_is1_core) THEN 572 center1 => exp_igr%core_centre(:, particle_set(p1)%shell_index) 573 ex1 => exp_igr%core_ex(:, particle_set(p1)%shell_index) 574 ey1 => exp_igr%core_ey(:, particle_set(p1)%shell_index) 575 ez1 => exp_igr%core_ez(:, particle_set(p1)%shell_index) 576 ELSE 577 center1 => exp_igr%centre(:, p1) 578 ex1 => exp_igr%ex(:, p1) 579 ey1 => exp_igr%ey(:, p1) 580 ez1 => exp_igr%ez(:, p1) 581 END IF 582 583 CALL dg_get_strucfac(box%hmat, r1, grid_s%npts, grid_b%npts, center1, & 584 exp_igr%lb, ex1, ey1, ez1) 585 586 IF (p2 /= 0) THEN 587 IF (my_is2_core) THEN 588 r2 = core_particle_set(particle_set(p2)%shell_index)%r 589 ELSE 590 r2 = particle_set(p2)%r 591 END IF 592 atomic_kind => particle_set(p2)%atomic_kind 593 IF (my_is2_core) THEN 594 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell) 595 q2 = shell%charge_core 596 ELSE IF (my_is2_shell) THEN 597 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell) 598 q2 = shell%charge_shell 599 ELSE 600 CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q2) 601 END IF 602 IF (use_charge_array) q2 = charges(p2) 603 604 IF (my_is2_shell) THEN 605 center2 => exp_igr%shell_centre(:, p2) 606 ex2 => exp_igr%shell_ex(:, p2) 607 ey2 => exp_igr%shell_ey(:, p2) 608 ez2 => exp_igr%shell_ez(:, p2) 609 ELSEIF (my_is2_core) THEN 610 center2 => exp_igr%core_centre(:, particle_set(p2)%shell_index) 611 ex2 => exp_igr%core_ex(:, particle_set(p2)%shell_index) 612 ey2 => exp_igr%core_ey(:, particle_set(p2)%shell_index) 613 ez2 => exp_igr%core_ez(:, particle_set(p2)%shell_index) 614 ELSE 615 center2 => exp_igr%centre(:, p2) 616 ex2 => exp_igr%ex(:, p2) 617 ey2 => exp_igr%ey(:, p2) 618 ez2 => exp_igr%ez(:, p2) 619 END IF 620 CALL dg_get_strucfac(box%hmat, r2, grid_s%npts, grid_b%npts, center2, & 621 exp_igr%lb, ex2, ey2, ez2) 622 END IF 623 624 IF (p2 == 0) THEN 625 CALL dg_get_patch(rho0, rhos1, q1, ex1, ey1, ez1) 626 ELSE 627 CALL dg_get_patch(rho0, rhos1, rhos2, q1, q2, ex1, ey1, ez1, ex2, ey2, ez2) 628 END IF 629 630 END SUBROUTINE get_patch 631 632! ************************************************************************************************** 633!> \brief ... 634!> \param dg ... 635!> \param particle_set ... 636!> \param exp_igr ... 637!> \param p1 ... 638!> \param p2 ... 639!> \param rhos1 ... 640!> \param rhos2 ... 641!> \param is1_core ... 642!> \param is2_core ... 643!> \param is1_shell ... 644!> \param is2_shell ... 645!> \param charges ... 646! ************************************************************************************************** 647 SUBROUTINE get_patch_again(dg, particle_set, exp_igr, p1, p2, rhos1, rhos2, is1_core, & 648 is2_core, is1_shell, is2_shell, charges) 649 650 TYPE(dg_type), POINTER :: dg 651 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 652 TYPE(structure_factor_type) :: exp_igr 653 INTEGER, INTENT(IN) :: p1, p2 654 TYPE(pw_p_type) :: rhos1, rhos2 655 LOGICAL, OPTIONAL :: is1_core, is2_core, is1_shell, is2_shell 656 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges 657 658 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_patch_again', & 659 routineP = moduleN//':'//routineN 660 661 COMPLEX(KIND=dp), DIMENSION(:), POINTER :: ex1, ex2, ey1, ey2, ez1, ez2 662 LOGICAL :: my_is1_core, my_is1_shell, my_is2_core, & 663 my_is2_shell, use_charge_array 664 REAL(KIND=dp) :: q1, q2 665 TYPE(atomic_kind_type), POINTER :: atomic_kind 666 TYPE(dg_rho0_type), POINTER :: dg_rho0 667 TYPE(pw_p_type), POINTER :: rho0 668 TYPE(shell_kind_type), POINTER :: shell 669 670 NULLIFY (shell) 671 use_charge_array = PRESENT(charges) 672 IF (use_charge_array) use_charge_array = ASSOCIATED(charges) 673 my_is1_core = .FALSE. 674 my_is2_core = .FALSE. 675 IF (PRESENT(is1_core)) my_is1_core = is1_core 676 IF (PRESENT(is2_core)) my_is2_core = is2_core 677 my_is1_shell = .FALSE. 678 my_is2_shell = .FALSE. 679 IF (PRESENT(is1_shell)) my_is1_shell = is1_shell 680 IF (PRESENT(is2_shell)) my_is2_shell = is2_shell 681 682 CALL dg_get(dg, dg_rho0=dg_rho0) 683 rho0 => dg_rho0%density 684 685 atomic_kind => particle_set(p1)%atomic_kind 686 IF (my_is1_core) THEN 687 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell) 688 q1 = shell%charge_core 689 ELSE IF (my_is1_shell) THEN 690 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell) 691 q1 = shell%charge_shell 692 ELSE 693 CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q1) 694 END IF 695 IF (use_charge_array) q1 = charges(p1) 696 IF (my_is1_core) THEN 697 ex1 => exp_igr%core_ex(:, particle_set(p1)%shell_index) 698 ey1 => exp_igr%core_ey(:, particle_set(p1)%shell_index) 699 ez1 => exp_igr%core_ez(:, particle_set(p1)%shell_index) 700 ELSEIF (my_is1_shell) THEN 701 ex1 => exp_igr%shell_ex(:, p1) 702 ey1 => exp_igr%shell_ey(:, p1) 703 ez1 => exp_igr%shell_ez(:, p1) 704 ELSE 705 ex1 => exp_igr%ex(:, p1) 706 ey1 => exp_igr%ey(:, p1) 707 ez1 => exp_igr%ez(:, p1) 708 END IF 709 710 IF (p2 /= 0) THEN 711 atomic_kind => particle_set(p2)%atomic_kind 712 IF (my_is2_core) THEN 713 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell) 714 q2 = shell%charge_core 715 ELSE IF (my_is2_shell) THEN 716 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell) 717 q2 = shell%charge_shell 718 ELSE 719 CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q2) 720 END IF 721 IF (use_charge_array) q2 = charges(p2) 722 IF (my_is2_core) THEN 723 ex2 => exp_igr%core_ex(:, particle_set(p2)%shell_index) 724 ey2 => exp_igr%core_ey(:, particle_set(p2)%shell_index) 725 ez2 => exp_igr%core_ez(:, particle_set(p2)%shell_index) 726 ELSEIF (my_is2_shell) THEN 727 ex2 => exp_igr%shell_ex(:, p2) 728 ey2 => exp_igr%shell_ey(:, p2) 729 ez2 => exp_igr%shell_ez(:, p2) 730 ELSE 731 ex2 => exp_igr%ex(:, p2) 732 ey2 => exp_igr%ey(:, p2) 733 ez2 => exp_igr%ez(:, p2) 734 END IF 735 END IF 736 737 IF (p2 == 0) THEN 738 CALL dg_get_patch(rho0, rhos1, q1, ex1, ey1, ez1) 739 ELSE 740 CALL dg_get_patch(rho0, rhos1, rhos2, q1, q2, & 741 ex1, ey1, ez1, ex2, ey2, ez2) 742 END IF 743 744 END SUBROUTINE get_patch_again 745 746END MODULE pme 747 748