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