(git:e7e05ae)
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 ! **************************************************************************************************
13 MODULE 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 
70 CONTAINS
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 
737 END MODULE pme
738 
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
dg_types::dg_type
Definition: dg_types.F:23