(git:ed6f26b)
Loading...
Searching...
No Matches
pme.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \par History
10!> JGH (21-Mar-2001) : Complete rewrite
11!> \author CJM and APSI
12! **************************************************************************************************
13MODULE pme
14
15 USE atomic_kind_types, ONLY: atomic_kind_type,&
16 get_atomic_kind
17 USE atprop_types, ONLY: atprop_type
18 USE bibliography, ONLY: cite_reference,&
19 darden1993
20 USE cell_types, ONLY: cell_type
21 USE dg_rho0_types, ONLY: dg_rho0_type
22 USE dg_types, ONLY: dg_get,&
23 dg_type
24 USE dgs, ONLY: dg_get_patch,&
25 dg_get_strucfac,&
26 dg_sum_patch,&
27 dg_sum_patch_force_1d,&
28 dg_sum_patch_force_3d
29 USE ewald_environment_types, ONLY: ewald_env_get,&
30 ewald_environment_type
31 USE ewald_pw_types, ONLY: ewald_pw_get,&
32 ewald_pw_type
33 USE kinds, ONLY: dp
34 USE mathconstants, ONLY: fourpi
35 USE message_passing, ONLY: mp_comm_type
36 USE particle_types, ONLY: particle_type
37 USE pme_tools, ONLY: get_center,&
38 set_list
39 USE pw_grid_types, ONLY: pw_grid_type
40 USE pw_methods, ONLY: pw_integral_a2b,&
41 pw_transfer
42 USE pw_poisson_methods, ONLY: pw_poisson_solve
43 USE pw_poisson_types, ONLY: pw_poisson_type
44 USE pw_pool_types, ONLY: pw_pool_type
45 USE pw_types, ONLY: pw_c1d_gs_type,&
46 pw_r3d_rs_type
47 USE realspace_grid_types, ONLY: realspace_grid_desc_type,&
48 realspace_grid_type,&
49 rs_grid_create,&
50 rs_grid_release,&
51 rs_grid_set_box,&
52 rs_grid_zero,&
53 transfer_pw2rs,&
54 transfer_rs2pw
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'
112
113 INTEGER :: handle, i, ipart, j, npart, nshell, p1, &
114 p2
115 LOGICAL :: is1_core, is2_core
116 REAL(kind=dp) :: alpha, dvols, fat1, ffa
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(mp_comm_type) :: group
122 TYPE(pw_c1d_gs_type), DIMENSION(3) :: dphi_g
123 TYPE(pw_grid_type), POINTER :: grid_b, grid_s
124 TYPE(pw_poisson_type), POINTER :: poisson_env
125 TYPE(pw_pool_type), POINTER :: pw_big_pool, pw_small_pool
126 TYPE(pw_r3d_rs_type) :: phi_r, rhob_r, rhos1, rhos2
127 TYPE(realspace_grid_desc_type), POINTER :: rs_desc
128 TYPE(realspace_grid_type), DIMENSION(3) :: 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)
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_small_pool%create_pw(rhos1)
163 CALL pw_small_pool%create_pw(rhos2)
164
165 ALLOCATE (rden)
166 CALL rs_grid_create(rden, rs_desc)
167 CALL rs_grid_set_box(grid_b, rs=rden)
168 CALL rs_grid_zero(rden)
169
170 cpassert(ASSOCIATED(box))
171
172 IF (rden%desc%parallel .AND. rden%desc%distributed) THEN
173 CALL get_center(particle_set, box, exp_igr%centre, exp_igr%delta, grid_b%npts, 1)
174 END IF
175 IF (PRESENT(shell_particle_set) .AND. rden%desc%parallel .AND. rden%desc%distributed) THEN
176 CALL get_center(shell_particle_set, box, exp_igr%shell_centre, exp_igr%shell_delta, grid_b%npts, 1)
177 CALL get_center(core_particle_set, box, exp_igr%core_centre, exp_igr%core_delta, grid_b%npts, 1)
178 END IF
179
180 !-------------- DENSITY CALCULATION ----------------
181
182 ipart = 0
183 DO
184
185 CALL set_list(particle_set, npart, exp_igr%centre, p1, rden, ipart, exp_igr%core_centre)
186 CALL set_list(particle_set, npart, exp_igr%centre, p2, rden, ipart, exp_igr%core_centre)
187 IF (p1 == 0 .AND. p2 == 0) EXIT
188
189 is1_core = (particle_set(p1)%shell_index /= 0)
190 IF (p2 /= 0) THEN
191 is2_core = (particle_set(p2)%shell_index /= 0)
192 ELSE
193 is2_core = .false.
194 END IF
195
196 ! calculate function on small boxes (we use double packing in FFT)
197 IF (is1_core .OR. is2_core) THEN
198 CALL get_patch(dg, particle_set, exp_igr, box, p1, p2, grid_b, grid_s, &
199 rhos1, rhos2, is1_core=is1_core, is2_core=is2_core, &
200 core_particle_set=core_particle_set, charges=charges)
201
202 ! add boxes to real space grid (big box)
203 IF (is1_core) THEN
204 CALL dg_sum_patch(rden, rhos1, exp_igr%core_centre(:, particle_set(p1)%shell_index))
205 ELSE
206 CALL dg_sum_patch(rden, rhos1, exp_igr%centre(:, p1))
207 END IF
208 IF (p2 /= 0 .AND. is2_core) THEN
209 CALL dg_sum_patch(rden, rhos2, exp_igr%core_centre(:, particle_set(p2)%shell_index))
210 ELSE IF (p2 /= 0) THEN
211 CALL dg_sum_patch(rden, rhos2, exp_igr%centre(:, p2))
212 END IF
213 ELSE
214 CALL get_patch(dg, particle_set, exp_igr, box, p1, p2, grid_b, grid_s, &
215 rhos1, rhos2, charges=charges)
216 ! add boxes to real space grid (big box)
217 CALL dg_sum_patch(rden, rhos1, exp_igr%centre(:, p1))
218 IF (p2 /= 0) CALL dg_sum_patch(rden, rhos2, exp_igr%centre(:, p2))
219 END IF
220
221 END DO
222 IF (PRESENT(shell_particle_set)) THEN
223 ipart = 0
224 DO
225 CALL set_list(shell_particle_set, nshell, exp_igr%shell_centre, p1, rpot, ipart)
226 CALL set_list(shell_particle_set, nshell, exp_igr%shell_centre, p2, rpot, ipart)
227 IF (p1 == 0 .AND. p2 == 0) EXIT
228 ! calculate function on small boxes (we use double packing in FFT)
229 CALL get_patch(dg, shell_particle_set, exp_igr, box, p1, p2, grid_b, grid_s, &
230 rhos1, rhos2, is1_shell=.true., is2_shell=.true., charges=charges)
231 ! add boxes to real space grid (big box)
232 CALL dg_sum_patch(rpot, rhos1, exp_igr%shell_centre(:, p1))
233 IF (p2 /= 0) CALL dg_sum_patch(rpot, rhos2, exp_igr%shell_centre(:, p2))
234 END DO
235 END IF
236
237 CALL pw_big_pool%create_pw(rhob_r)
238 CALL transfer_rs2pw(rden, rhob_r)
239
240 !-------------- ELECTROSTATIC CALCULATION -----------
241
242 ! allocate intermediate arrays
243 DO i = 1, 3
244 CALL pw_big_pool%create_pw(dphi_g(i))
245 END DO
246 CALL pw_big_pool%create_pw(phi_r)
247
248 CALL pw_poisson_solve(poisson_env, rhob_r, vg_coulomb, phi_r, dphi_g, h_stress)
249
250 ! atomic energies
251 IF (atprop%energy) THEN
252 dvols = rhos1%pw_grid%dvol
253 ALLOCATE (rpot)
254 CALL rs_grid_create(rpot, rs_desc)
255 CALL transfer_pw2rs(rpot, phi_r)
256 ipart = 0
257 DO
258 CALL set_list(particle_set, npart, exp_igr%centre, p1, rden, ipart, exp_igr%core_centre)
259 CALL set_list(particle_set, npart, exp_igr%centre, p2, rden, ipart, exp_igr%core_centre)
260 IF (p1 == 0 .AND. p2 == 0) EXIT
261 ! integrate box and potential
262 CALL get_patch(dg, particle_set, exp_igr, box, p1, p2, grid_b, grid_s, &
263 rhos1, rhos2, charges=charges)
264 ! add boxes to real space grid (big box)
265 CALL dg_sum_patch_force_1d(rpot, rhos1, exp_igr%centre(:, p1), fat1)
266 IF (atprop%energy) THEN
267 atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols
268 END IF
269 IF (p2 /= 0) THEN
270 CALL dg_sum_patch_force_1d(rpot, rhos2, exp_igr%centre(:, p2), fat1)
271 IF (atprop%energy) THEN
272 atprop%atener(p2) = atprop%atener(p2) + 0.5_dp*fat1*dvols
273 END IF
274 END IF
275 END DO
276 CALL rs_grid_release(rpot)
277 DEALLOCATE (rpot)
278 END IF
279
280 CALL pw_big_pool%give_back_pw(phi_r)
281
282 !---------- END OF ELECTROSTATIC CALCULATION --------
283
284 !------------- STRESS TENSOR CALCULATION ------------
285
286 IF ((use_virial) .AND. (PRESENT(pv_g))) THEN
287 DO i = 1, 3
288 DO j = i, 3
289 f_stress(i, j) = pw_integral_a2b(dphi_g(i), dphi_g(j))
290 f_stress(j, i) = f_stress(i, j)
291 END DO
292 END DO
293 ffa = (1.0_dp/fourpi)*(0.5_dp/dg_rho0%zet(1))**2
294 f_stress = -ffa*f_stress
295 pv_g = h_stress + f_stress
296 END IF
297
298 !--------END OF STRESS TENSOR CALCULATION -----------
299
300 DO i = 1, 3
301 CALL rs_grid_create(drpot(i), rs_desc)
302 CALL rs_grid_set_box(grid_b, rs=drpot(i))
303 CALL pw_transfer(dphi_g(i), rhob_r)
304 CALL pw_big_pool%give_back_pw(dphi_g(i))
305 CALL transfer_pw2rs(drpot(i), rhob_r)
306 END DO
307
308 CALL pw_big_pool%give_back_pw(rhob_r)
309
310 !----------------- FORCE CALCULATION ----------------
311
312 ! initialize the forces
313 IF (PRESENT(fg_coulomb)) THEN
314 fg_coulomb = 0.0_dp
315 dvols = rhos1%pw_grid%dvol
316
317 ipart = 0
318 DO
319
320 CALL set_list(particle_set, npart, exp_igr%centre, p1, rden, ipart, exp_igr%core_centre)
321 CALL set_list(particle_set, npart, exp_igr%centre, p2, rden, ipart, exp_igr%core_centre)
322 IF (p1 == 0 .AND. p2 == 0) EXIT
323
324 is1_core = (particle_set(p1)%shell_index /= 0)
325 IF (p2 /= 0) THEN
326 is2_core = (particle_set(p2)%shell_index /= 0)
327 ELSE
328 is2_core = .false.
329 END IF
330
331 ! calculate function on small boxes (we use double packing in FFT)
332
333 CALL get_patch_again(dg, particle_set, exp_igr, p1, p2, rhos1, rhos2, &
334 is1_core=is1_core, is2_core=is2_core, charges=charges)
335
336 ! sum boxes on real space grids (big box)
337 IF (is1_core) THEN
338 CALL dg_sum_patch_force_3d(drpot, rhos1, &
339 exp_igr%core_centre(:, particle_set(p1)%shell_index), fat)
340 fgcore_coulomb(1, particle_set(p1)%shell_index) = &
341 fgcore_coulomb(1, particle_set(p1)%shell_index) - fat(1)*dvols
342 fgcore_coulomb(2, particle_set(p1)%shell_index) = &
343 fgcore_coulomb(2, particle_set(p1)%shell_index) - fat(2)*dvols
344 fgcore_coulomb(3, particle_set(p1)%shell_index) = &
345 fgcore_coulomb(3, particle_set(p1)%shell_index) - fat(3)*dvols
346 ELSE
347 CALL dg_sum_patch_force_3d(drpot, rhos1, exp_igr%centre(:, p1), fat)
348 fg_coulomb(1, p1) = fg_coulomb(1, p1) - fat(1)*dvols
349 fg_coulomb(2, p1) = fg_coulomb(2, p1) - fat(2)*dvols
350 fg_coulomb(3, p1) = fg_coulomb(3, p1) - fat(3)*dvols
351 END IF
352 IF (p2 /= 0 .AND. is2_core) THEN
353 CALL dg_sum_patch_force_3d(drpot, rhos1, &
354 exp_igr%core_centre(:, particle_set(p2)%shell_index), fat)
355 fgcore_coulomb(1, particle_set(p2)%shell_index) = &
356 fgcore_coulomb(1, particle_set(p2)%shell_index) - fat(1)*dvols
357 fgcore_coulomb(2, particle_set(p2)%shell_index) = &
358 fgcore_coulomb(2, particle_set(p2)%shell_index) - fat(2)*dvols
359 fgcore_coulomb(3, particle_set(p2)%shell_index) = &
360 fgcore_coulomb(3, particle_set(p2)%shell_index) - fat(3)*dvols
361 ELSEIF (p2 /= 0) THEN
362 CALL dg_sum_patch_force_3d(drpot, rhos2, exp_igr%centre(:, p2), fat)
363 fg_coulomb(1, p2) = fg_coulomb(1, p2) - fat(1)*dvols
364 fg_coulomb(2, p2) = fg_coulomb(2, p2) - fat(2)*dvols
365 fg_coulomb(3, p2) = fg_coulomb(3, p2) - fat(3)*dvols
366 END IF
367
368 END DO
369 END IF
370 IF (PRESENT(fgshell_coulomb)) THEN
371 fgshell_coulomb = 0.0_dp
372 dvols = rhos1%pw_grid%dvol
373
374 ipart = 0
375 DO
376 CALL set_list(shell_particle_set, nshell, exp_igr%shell_centre, p1, rden, ipart)
377 CALL set_list(shell_particle_set, nshell, exp_igr%shell_centre, p2, rden, ipart)
378 IF (p1 == 0 .AND. p2 == 0) EXIT
379
380 ! calculate function on small boxes (we use double packing in FFT)
381 CALL get_patch_again(dg, shell_particle_set, exp_igr, p1, p2, rhos1, rhos2, &
382 is1_shell=.true., is2_shell=.true., charges=charges)
383
384 ! sum boxes on real space grids (big box)
385 CALL dg_sum_patch_force_3d(drpot, rhos1, exp_igr%shell_centre(:, p1), fat)
386 fgshell_coulomb(1, p1) = fgshell_coulomb(1, p1) - fat(1)*dvols
387 fgshell_coulomb(2, p1) = fgshell_coulomb(2, p1) - fat(2)*dvols
388 fgshell_coulomb(3, p1) = fgshell_coulomb(3, p1) - fat(3)*dvols
389 IF (p2 /= 0) THEN
390 CALL dg_sum_patch_force_3d(drpot, rhos2, exp_igr%shell_centre(:, p2), fat)
391 fgshell_coulomb(1, p2) = fgshell_coulomb(1, p2) - fat(1)*dvols
392 fgshell_coulomb(2, p2) = fgshell_coulomb(2, p2) - fat(2)*dvols
393 fgshell_coulomb(3, p2) = fgshell_coulomb(3, p2) - fat(3)*dvols
394 END IF
395 END DO
396
397 END IF
398 !--------------END OF FORCE CALCULATION -------------
399
400 !------------------CLEANING UP ----------------------
401
402 CALL rs_grid_release(rden)
403 DEALLOCATE (rden)
404 DO i = 1, 3
405 CALL rs_grid_release(drpot(i))
406 END DO
407
408 CALL pw_small_pool%give_back_pw(rhos1)
409 CALL pw_small_pool%give_back_pw(rhos2)
410 CALL structure_factor_deallocate(exp_igr)
411
412 CALL timestop(handle)
413
414 END SUBROUTINE pme_evaluate
415
416! **************************************************************************************************
417!> \brief Calculates local density in a small box
418!> \param dg ...
419!> \param particle_set ...
420!> \param exp_igr ...
421!> \param box ...
422!> \param p1 ...
423!> \param p2 ...
424!> \param grid_b ...
425!> \param grid_s ...
426!> \param rhos1 ...
427!> \param rhos2 ...
428!> \param is1_core ...
429!> \param is2_core ...
430!> \param is1_shell ...
431!> \param is2_shell ...
432!> \param core_particle_set ...
433!> \param charges ...
434!> \par History
435!> JGH (23-Mar-2001) : Switch to integer from particle list pointers
436!> \author JGH (21-Mar-2001)
437! **************************************************************************************************
438 SUBROUTINE get_patch(dg, particle_set, exp_igr, box, p1, p2, &
439 grid_b, grid_s, rhos1, rhos2, is1_core, is2_core, is1_shell, &
440 is2_shell, core_particle_set, charges)
441
442 TYPE(dg_type), POINTER :: dg
443 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
444 TYPE(structure_factor_type) :: exp_igr
445 TYPE(cell_type), POINTER :: box
446 INTEGER, INTENT(IN) :: p1, p2
447 TYPE(pw_grid_type), INTENT(IN) :: grid_b, grid_s
448 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rhos1, rhos2
449 LOGICAL, OPTIONAL :: is1_core, is2_core, is1_shell, is2_shell
450 TYPE(particle_type), DIMENSION(:), OPTIONAL, &
451 POINTER :: core_particle_set
452 REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: charges
453
454 COMPLEX(KIND=dp), DIMENSION(:), POINTER :: ex1, ex2, ey1, ey2, ez1, ez2
455 INTEGER, DIMENSION(:), POINTER :: center1, center2
456 LOGICAL :: my_is1_core, my_is1_shell, my_is2_core, &
457 my_is2_shell, use_charge_array
458 REAL(kind=dp) :: q1, q2
459 REAL(kind=dp), DIMENSION(3) :: r1, r2
460 TYPE(atomic_kind_type), POINTER :: atomic_kind
461 TYPE(dg_rho0_type), POINTER :: dg_rho0
462 TYPE(pw_r3d_rs_type), POINTER :: rho0
463 TYPE(shell_kind_type), POINTER :: shell
464
465 NULLIFY (shell)
466 use_charge_array = .false.
467 IF (PRESENT(charges)) use_charge_array = ASSOCIATED(charges)
468 my_is1_core = .false.
469 my_is2_core = .false.
470 IF (PRESENT(is1_core)) my_is1_core = is1_core
471 IF (PRESENT(is2_core)) my_is2_core = is2_core
472 IF (my_is1_core .OR. my_is2_core) THEN
473 cpassert(PRESENT(core_particle_set))
474 END IF
475 my_is1_shell = .false.
476 my_is2_shell = .false.
477 IF (PRESENT(is1_shell)) my_is1_shell = is1_shell
478 IF (PRESENT(is2_shell)) my_is2_shell = is2_shell
479 IF (my_is1_core .AND. my_is1_shell) THEN
480 cpabort("Shell-model: cannot be core and shell simultaneously")
481 END IF
482
483 CALL dg_get(dg, dg_rho0=dg_rho0)
484 rho0 => dg_rho0%density
485
486 IF (my_is1_core) THEN
487 r1 = core_particle_set(particle_set(p1)%shell_index)%r
488 ELSE
489 r1 = particle_set(p1)%r
490 END IF
491 atomic_kind => particle_set(p1)%atomic_kind
492 IF (my_is1_core) THEN
493 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
494 q1 = shell%charge_core
495 ELSE IF (my_is1_shell) THEN
496 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
497 q1 = shell%charge_shell
498 ELSE
499 CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q1)
500 END IF
501 IF (use_charge_array) q1 = charges(p1)
502
503 IF (my_is1_shell) THEN
504 center1 => exp_igr%shell_centre(:, p1)
505 ex1 => exp_igr%shell_ex(:, p1)
506 ey1 => exp_igr%shell_ey(:, p1)
507 ez1 => exp_igr%shell_ez(:, p1)
508 ELSEIF (my_is1_core) THEN
509 center1 => exp_igr%core_centre(:, particle_set(p1)%shell_index)
510 ex1 => exp_igr%core_ex(:, particle_set(p1)%shell_index)
511 ey1 => exp_igr%core_ey(:, particle_set(p1)%shell_index)
512 ez1 => exp_igr%core_ez(:, particle_set(p1)%shell_index)
513 ELSE
514 center1 => exp_igr%centre(:, p1)
515 ex1 => exp_igr%ex(:, p1)
516 ey1 => exp_igr%ey(:, p1)
517 ez1 => exp_igr%ez(:, p1)
518 END IF
519
520 cpassert(ASSOCIATED(box))
521
522 CALL dg_get_strucfac(box%hmat, r1, grid_s%npts, grid_b%npts, center1, &
523 exp_igr%lb, ex1, ey1, ez1)
524
525 IF (p2 /= 0) THEN
526 IF (my_is2_core) THEN
527 r2 = core_particle_set(particle_set(p2)%shell_index)%r
528 ELSE
529 r2 = particle_set(p2)%r
530 END IF
531 atomic_kind => particle_set(p2)%atomic_kind
532 IF (my_is2_core) THEN
533 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
534 q2 = shell%charge_core
535 ELSE IF (my_is2_shell) THEN
536 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
537 q2 = shell%charge_shell
538 ELSE
539 CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q2)
540 END IF
541 IF (use_charge_array) q2 = charges(p2)
542
543 IF (my_is2_shell) THEN
544 center2 => exp_igr%shell_centre(:, p2)
545 ex2 => exp_igr%shell_ex(:, p2)
546 ey2 => exp_igr%shell_ey(:, p2)
547 ez2 => exp_igr%shell_ez(:, p2)
548 ELSEIF (my_is2_core) THEN
549 center2 => exp_igr%core_centre(:, particle_set(p2)%shell_index)
550 ex2 => exp_igr%core_ex(:, particle_set(p2)%shell_index)
551 ey2 => exp_igr%core_ey(:, particle_set(p2)%shell_index)
552 ez2 => exp_igr%core_ez(:, particle_set(p2)%shell_index)
553 ELSE
554 center2 => exp_igr%centre(:, p2)
555 ex2 => exp_igr%ex(:, p2)
556 ey2 => exp_igr%ey(:, p2)
557 ez2 => exp_igr%ez(:, p2)
558 END IF
559 CALL dg_get_strucfac(box%hmat, r2, grid_s%npts, grid_b%npts, center2, &
560 exp_igr%lb, ex2, ey2, ez2)
561 END IF
562
563 IF (p2 == 0) THEN
564 CALL dg_get_patch(rho0, rhos1, q1, ex1, ey1, ez1)
565 ELSE
566 CALL dg_get_patch(rho0, rhos1, rhos2, q1, q2, ex1, ey1, ez1, ex2, ey2, ez2)
567 END IF
568
569 END SUBROUTINE get_patch
570
571! **************************************************************************************************
572!> \brief ...
573!> \param dg ...
574!> \param particle_set ...
575!> \param exp_igr ...
576!> \param p1 ...
577!> \param p2 ...
578!> \param rhos1 ...
579!> \param rhos2 ...
580!> \param is1_core ...
581!> \param is2_core ...
582!> \param is1_shell ...
583!> \param is2_shell ...
584!> \param charges ...
585! **************************************************************************************************
586 SUBROUTINE get_patch_again(dg, particle_set, exp_igr, p1, p2, rhos1, rhos2, is1_core, &
587 is2_core, is1_shell, is2_shell, charges)
588
589 TYPE(dg_type), POINTER :: dg
590 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
591 TYPE(structure_factor_type) :: exp_igr
592 INTEGER, INTENT(IN) :: p1, p2
593 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rhos1, rhos2
594 LOGICAL, OPTIONAL :: is1_core, is2_core, is1_shell, is2_shell
595 REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: charges
596
597 COMPLEX(KIND=dp), DIMENSION(:), POINTER :: ex1, ex2, ey1, ey2, ez1, ez2
598 LOGICAL :: my_is1_core, my_is1_shell, my_is2_core, &
599 my_is2_shell, use_charge_array
600 REAL(kind=dp) :: q1, q2
601 TYPE(atomic_kind_type), POINTER :: atomic_kind
602 TYPE(dg_rho0_type), POINTER :: dg_rho0
603 TYPE(pw_r3d_rs_type), POINTER :: rho0
604 TYPE(shell_kind_type), POINTER :: shell
605
606 NULLIFY (shell)
607 use_charge_array = .false.
608 IF (PRESENT(charges)) use_charge_array = ASSOCIATED(charges)
609 my_is1_core = .false.
610 my_is2_core = .false.
611 IF (PRESENT(is1_core)) my_is1_core = is1_core
612 IF (PRESENT(is2_core)) my_is2_core = is2_core
613 my_is1_shell = .false.
614 my_is2_shell = .false.
615 IF (PRESENT(is1_shell)) my_is1_shell = is1_shell
616 IF (PRESENT(is2_shell)) my_is2_shell = is2_shell
617
618 CALL dg_get(dg, dg_rho0=dg_rho0)
619 rho0 => dg_rho0%density
620
621 atomic_kind => particle_set(p1)%atomic_kind
622 IF (my_is1_core) THEN
623 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
624 q1 = shell%charge_core
625 ELSE IF (my_is1_shell) THEN
626 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
627 q1 = shell%charge_shell
628 ELSE
629 CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q1)
630 END IF
631 IF (use_charge_array) q1 = charges(p1)
632 IF (my_is1_core) THEN
633 ex1 => exp_igr%core_ex(:, particle_set(p1)%shell_index)
634 ey1 => exp_igr%core_ey(:, particle_set(p1)%shell_index)
635 ez1 => exp_igr%core_ez(:, particle_set(p1)%shell_index)
636 ELSEIF (my_is1_shell) THEN
637 ex1 => exp_igr%shell_ex(:, p1)
638 ey1 => exp_igr%shell_ey(:, p1)
639 ez1 => exp_igr%shell_ez(:, p1)
640 ELSE
641 ex1 => exp_igr%ex(:, p1)
642 ey1 => exp_igr%ey(:, p1)
643 ez1 => exp_igr%ez(:, p1)
644 END IF
645
646 IF (p2 /= 0) THEN
647 atomic_kind => particle_set(p2)%atomic_kind
648 IF (my_is2_core) THEN
649 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
650 q2 = shell%charge_core
651 ELSE IF (my_is2_shell) THEN
652 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
653 q2 = shell%charge_shell
654 ELSE
655 CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q2)
656 END IF
657 IF (use_charge_array) q2 = charges(p2)
658 IF (my_is2_core) THEN
659 ex2 => exp_igr%core_ex(:, particle_set(p2)%shell_index)
660 ey2 => exp_igr%core_ey(:, particle_set(p2)%shell_index)
661 ez2 => exp_igr%core_ez(:, particle_set(p2)%shell_index)
662 ELSEIF (my_is2_shell) THEN
663 ex2 => exp_igr%shell_ex(:, p2)
664 ey2 => exp_igr%shell_ey(:, p2)
665 ez2 => exp_igr%shell_ez(:, p2)
666 ELSE
667 ex2 => exp_igr%ex(:, p2)
668 ey2 => exp_igr%ey(:, p2)
669 ez2 => exp_igr%ez(:, p2)
670 END IF
671 END IF
672
673 IF (p2 == 0) THEN
674 CALL dg_get_patch(rho0, rhos1, q1, ex1, ey1, ez1)
675 ELSE
676 CALL dg_get_patch(rho0, rhos1, rhos2, q1, q2, &
677 ex1, ey1, ez1, ex2, ey2, ez2)
678 END IF
679
680 END SUBROUTINE get_patch_again
681
682END MODULE pme
683
dgs::dg_get_patch
Definition dgs.F:64
dgs::dg_sum_patch_force_1d
Definition dgs.F:60
dgs::dg_sum_patch_force_3d
Definition dgs.F:56
dgs::dg_sum_patch
Definition dgs.F:52
pw_methods::pw_integral_a2b
Definition pw_methods.F:240
pw_methods::pw_transfer
Definition pw_methods.F:303
pw_poisson_methods::pw_poisson_solve
Definition pw_poisson_methods.F:78
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
atprop_types
Holds information on atomic properties.
Definition atprop_types.F:14
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::darden1993
integer, save, public darden1993
Definition bibliography.F:36
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
dg_rho0_types
Definition dg_rho0_types.F:12
dg_types
Definition dg_types.F:12
dg_types::dg_get
subroutine, public dg_get(dg, dg_rho0)
Get the dg_type.
Definition dg_types.F:44
dgs
Definition dgs.F:13
dgs::dg_get_strucfac
subroutine, public dg_get_strucfac(cell_hmat, r, npts_s, npts_b, centre, lb, ex, ey, ez)
...
Definition dgs.F:333
ewald_environment_types
Definition ewald_environment_types.F:14
ewald_environment_types::ewald_env_get
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
Definition ewald_environment_types.F:124
ewald_pw_types
pw_types
Definition ewald_pw_types.F:12
ewald_pw_types::ewald_pw_get
subroutine, public ewald_pw_get(ewald_pw, pw_big_pool, pw_small_pool, rs_desc, poisson_env, dg)
get the ewald_pw environment to the correct program.
Definition ewald_pw_types.F:320
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::fourpi
real(kind=dp), parameter, public fourpi
Definition mathconstants.F:113
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pme_tools
Tools common both to PME and SPME.
Definition pme_tools.F:15
pme_tools::get_center
subroutine, public get_center(part, box, center, delta, npts, n)
...
Definition pme_tools.F:131
pme_tools::set_list
subroutine, public set_list(part, npart, center, p1, rs, ipart, core_center)
...
Definition pme_tools.F:46
pme
Definition pme.F:13
pme::pme_evaluate
subroutine, public pme_evaluate(ewald_env, ewald_pw, box, particle_set, vg_coulomb, fg_coulomb, pv_g, shell_particle_set, core_particle_set, fgshell_coulomb, fgcore_coulomb, use_virial, charges, atprop)
...
Definition pme.F:96
pw_grid_types
Definition pw_grid_types.F:13
pw_methods
Definition pw_methods.F:25
pw_poisson_methods
Definition pw_poisson_methods.F:13
pw_poisson_types
functions related to the poisson solver on regular grids
Definition pw_poisson_types.F:22
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
realspace_grid_types
Definition realspace_grid_types.F:18
realspace_grid_types::rs_grid_create
subroutine, public rs_grid_create(rs, desc)
...
Definition realspace_grid_types.F:478
realspace_grid_types::transfer_pw2rs
subroutine, public transfer_pw2rs(rs, pw)
...
Definition realspace_grid_types.F:694
realspace_grid_types::rs_grid_set_box
subroutine, public rs_grid_set_box(pw_grid, rs)
Set box matrix info for real space grid This is needed for variable cell simulations.
Definition realspace_grid_types.F:1980
realspace_grid_types::transfer_rs2pw
subroutine, public transfer_rs2pw(rs, pw)
...
Definition realspace_grid_types.F:652
realspace_grid_types::rs_grid_release
subroutine, public rs_grid_release(rs_grid)
releases the given rs grid (see doc/ReferenceCounting.html)
Definition realspace_grid_types.F:2012
realspace_grid_types::rs_grid_zero
subroutine, public rs_grid_zero(rs)
Initialize grid to zero.
Definition realspace_grid_types.F:1905
shell_potential_types
Definition shell_potential_types.F:11
structure_factor_types
Definition structure_factor_types.F:12
structure_factors
Definition structure_factors.F:12
structure_factors::structure_factor_deallocate
subroutine, public structure_factor_deallocate(exp_igr)
...
Definition structure_factors.F:50
structure_factors::structure_factor_allocate
subroutine, public structure_factor_allocate(bds, nparts, exp_igr, allocate_centre, allocate_shell_e, allocate_shell_centre, nshell)
...
Definition structure_factors.F:91
structure_factors::structure_factor_init
subroutine, public structure_factor_init(exp_igr)
...
Definition structure_factors.F:34
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
atprop_types::atprop_type
type for the atomic properties
Definition atprop_types.F:31
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
dg_rho0_types::dg_rho0_type
Type for Gaussian Densities type = type of gaussian (PME) grid = grid number gcc = Gaussian contracti...
Definition dg_rho0_types.F:39
dg_types::dg_type
Definition dg_types.F:23
ewald_environment_types::ewald_environment_type
to build arrays of pointers
Definition ewald_environment_types.F:56
ewald_pw_types::ewald_pw_type
Definition ewald_pw_types.F:63
message_passing::mp_comm_type
Definition message_passing.F:189
particle_types::particle_type
Definition particle_types.F:35
pw_grid_types::pw_grid_type
Definition pw_grid_types.F:52
pw_poisson_types::pw_poisson_type
environment for the poisson solver
Definition pw_poisson_types.F:122
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
realspace_grid_types::realspace_grid_desc_type
Definition realspace_grid_types.F:90
realspace_grid_types::realspace_grid_type
Definition realspace_grid_types.F:136
shell_potential_types::shell_kind_type
Define the shell type.
Definition shell_potential_types.F:28
structure_factor_types::structure_factor_type
Definition structure_factor_types.F:23