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