Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 : mp_para_env_type
32 : USE particle_types, ONLY: particle_type
33 : USE pme_tools, ONLY: get_center,&
34 : set_list
35 : USE pw_grid_types, ONLY: pw_grid_type
36 : USE pw_grids, ONLY: get_pw_grid_info
37 : USE pw_methods, ONLY: pw_integral_a2b,&
38 : pw_multiply_with,&
39 : pw_transfer
40 : USE pw_poisson_methods, ONLY: pw_poisson_rebuild,&
41 : pw_poisson_solve
42 : USE pw_poisson_types, ONLY: greens_fn_type,&
43 : pw_poisson_type
44 : USE pw_pool_types, ONLY: pw_pool_type
45 : USE pw_types, ONLY: pw_c1d_gs_type,&
46 : pw_r3d_rs_type
47 : USE realspace_grid_types, ONLY: realspace_grid_desc_type,&
48 : realspace_grid_type,&
49 : rs_grid_create,&
50 : rs_grid_release,&
51 : rs_grid_set_box,&
52 : rs_grid_zero,&
53 : transfer_pw2rs,&
54 : transfer_rs2pw
55 : USE shell_potential_types, ONLY: shell_kind_type
56 : #include "./base/base_uses.f90"
57 :
58 : IMPLICIT NONE
59 :
60 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'spme'
61 :
62 : PRIVATE
63 : PUBLIC :: spme_evaluate, spme_potential, spme_forces, spme_virial, get_patch
64 :
65 : INTERFACE get_patch
66 : MODULE PROCEDURE get_patch_a, get_patch_b
67 : END INTERFACE
68 :
69 : CONTAINS
70 :
71 : ! **************************************************************************************************
72 : !> \brief ...
73 : !> \param ewald_env ...
74 : !> \param ewald_pw ...
75 : !> \param box ...
76 : !> \param particle_set ...
77 : !> \param fg_coulomb ...
78 : !> \param vg_coulomb ...
79 : !> \param pv_g ...
80 : !> \param shell_particle_set ...
81 : !> \param core_particle_set ...
82 : !> \param fgshell_coulomb ...
83 : !> \param fgcore_coulomb ...
84 : !> \param use_virial ...
85 : !> \param charges ...
86 : !> \param atprop ...
87 : !> \par History
88 : !> JGH (03-May-2001) : SPME with charge definition
89 : !> \author JGH (21-Mar-2001)
90 : ! **************************************************************************************************
91 31928 : SUBROUTINE spme_evaluate(ewald_env, ewald_pw, box, particle_set, &
92 31928 : fg_coulomb, vg_coulomb, pv_g, shell_particle_set, core_particle_set, &
93 31928 : fgshell_coulomb, fgcore_coulomb, use_virial, charges, atprop)
94 :
95 : TYPE(ewald_environment_type), POINTER :: ewald_env
96 : TYPE(ewald_pw_type), POINTER :: ewald_pw
97 : TYPE(cell_type), POINTER :: box
98 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
99 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: fg_coulomb
100 : REAL(KIND=dp), INTENT(OUT) :: vg_coulomb
101 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: pv_g
102 : TYPE(particle_type), DIMENSION(:), OPTIONAL, &
103 : POINTER :: shell_particle_set, core_particle_set
104 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), &
105 : OPTIONAL :: fgshell_coulomb, fgcore_coulomb
106 : LOGICAL, INTENT(IN) :: use_virial
107 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges
108 : TYPE(atprop_type), POINTER :: atprop
109 :
110 : CHARACTER(len=*), PARAMETER :: routineN = 'spme_evaluate'
111 :
112 : INTEGER :: handle, i, ipart, j, n, ncore, npart, &
113 : nshell, o_spline, p1, p1_shell
114 63856 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: center, core_center, shell_center
115 : INTEGER, DIMENSION(3) :: npts
116 : LOGICAL :: do_shell
117 : REAL(KIND=dp) :: alpha, dvols, fat1, ffa
118 31928 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: core_delta, delta, shell_delta
119 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rhos
120 : REAL(KIND=dp), DIMENSION(3) :: fat
121 : REAL(KIND=dp), DIMENSION(3, 3) :: f_stress, h_stress
122 : TYPE(greens_fn_type), POINTER :: green
123 : TYPE(mp_comm_type) :: group
124 127712 : TYPE(pw_c1d_gs_type), DIMENSION(3) :: dphi_g
125 : TYPE(pw_c1d_gs_type), POINTER :: phi_g, rhob_g
126 : TYPE(pw_grid_type), POINTER :: grid_spme
127 : TYPE(pw_poisson_type), POINTER :: poisson_env
128 : TYPE(pw_pool_type), POINTER :: pw_pool
129 : TYPE(pw_r3d_rs_type), POINTER :: rhob_r
130 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
131 1149408 : TYPE(realspace_grid_type), DIMENSION(3) :: drpot
132 : TYPE(realspace_grid_type), POINTER :: rden, rpot
133 :
134 31928 : NULLIFY (grid_spme, green, poisson_env, phi_g, rhob_g, rhob_r, pw_pool, rden, rpot)
135 31928 : CALL timeset(routineN, handle)
136 31928 : CALL cite_reference(Essmann1995)
137 :
138 : !-------------- INITIALISATION ---------------------
139 : CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, &
140 31928 : group=group)
141 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, &
142 31928 : poisson_env=poisson_env)
143 31928 : CALL pw_poisson_rebuild(poisson_env)
144 31928 : green => poisson_env%green_fft
145 31928 : grid_spme => pw_pool%pw_grid
146 :
147 31928 : npart = SIZE(particle_set)
148 :
149 31928 : CALL get_pw_grid_info(grid_spme, npts=npts, dvol=dvols)
150 :
151 31928 : n = o_spline
152 159640 : ALLOCATE (rhos(n, n, n))
153 670488 : ALLOCATE (rden)
154 31928 : CALL rs_grid_create(rden, rs_desc)
155 31928 : CALL rs_grid_set_box(grid_spme, rs=rden)
156 31928 : CALL rs_grid_zero(rden)
157 :
158 159640 : ALLOCATE (center(3, npart), delta(3, npart))
159 31928 : CALL get_center(particle_set, box, center, delta, npts, n)
160 :
161 31928 : IF (PRESENT(shell_particle_set) .AND. (PRESENT(core_particle_set))) THEN
162 11732 : CPASSERT(ASSOCIATED(shell_particle_set))
163 11732 : CPASSERT(ASSOCIATED(core_particle_set))
164 11732 : nshell = SIZE(shell_particle_set)
165 11732 : ncore = SIZE(core_particle_set)
166 11732 : CPASSERT(nshell == ncore)
167 58660 : ALLOCATE (shell_center(3, nshell), shell_delta(3, nshell))
168 11732 : CALL get_center(shell_particle_set, box, shell_center, shell_delta, npts, n)
169 35196 : ALLOCATE (core_center(3, nshell), core_delta(3, nshell))
170 11732 : CALL get_center(core_particle_set, box, core_center, core_delta, npts, n)
171 : END IF
172 :
173 : !-------------- DENSITY CALCULATION ----------------
174 31928 : ipart = 0
175 : ! Particles
176 : DO
177 3618054 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
178 3618054 : IF (p1 == 0) EXIT
179 :
180 3586126 : do_shell = (particle_set(p1)%shell_index /= 0)
181 3586126 : IF (do_shell) CYCLE
182 : ! calculate function on small boxes
183 : CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
184 3156071 : is_shell=.FALSE., unit_charge=.FALSE., charges=charges)
185 :
186 : ! add boxes to real space grid (big box)
187 3618054 : CALL dg_sum_patch(rden, rhos, center(:, p1))
188 : END DO
189 : ! Shell-Model
190 31928 : IF (PRESENT(shell_particle_set) .AND. PRESENT(core_particle_set)) THEN
191 11732 : ipart = 0
192 441787 : DO
193 : ! calculate function on small boxes
194 : CALL set_list(shell_particle_set, nshell, shell_center, p1_shell, &
195 441787 : rden, ipart)
196 441787 : IF (p1_shell == 0) EXIT
197 : CALL get_patch(shell_particle_set, shell_delta, green, p1_shell, rhos, &
198 430055 : is_core=.FALSE., is_shell=.TRUE., unit_charge=.FALSE.)
199 :
200 : ! add boxes to real space grid (big box)
201 441787 : CALL dg_sum_patch(rden, rhos, shell_center(:, p1_shell))
202 : END DO
203 11732 : ipart = 0
204 441787 : DO
205 : ! calculate function on small boxes
206 : CALL set_list(core_particle_set, nshell, core_center, p1_shell, &
207 441787 : rden, ipart)
208 441787 : IF (p1_shell == 0) EXIT
209 : CALL get_patch(core_particle_set, core_delta, green, p1_shell, rhos, &
210 430055 : is_core=.TRUE., is_shell=.FALSE., unit_charge=.FALSE.)
211 :
212 : ! add boxes to real space grid (big box)
213 441787 : CALL dg_sum_patch(rden, rhos, core_center(:, p1_shell))
214 : END DO
215 : END IF
216 : !----------- END OF DENSITY CALCULATION -------------
217 :
218 31928 : NULLIFY (rhob_r)
219 31928 : ALLOCATE (rhob_r)
220 31928 : CALL pw_pool%create_pw(rhob_r)
221 31928 : CALL transfer_rs2pw(rden, rhob_r)
222 : ! transform density to G space and add charge function
223 31928 : NULLIFY (rhob_g)
224 31928 : ALLOCATE (rhob_g)
225 31928 : CALL pw_pool%create_pw(rhob_g)
226 31928 : CALL pw_transfer(rhob_r, rhob_g)
227 : ! update charge function
228 31928 : CALL pw_multiply_with(rhob_g, green%p3m_charge)
229 :
230 : !-------------- ELECTROSTATIC CALCULATION -----------
231 : ! allocate intermediate arrays
232 127712 : DO i = 1, 3
233 127712 : CALL pw_pool%create_pw(dphi_g(i))
234 : END DO
235 31928 : NULLIFY (phi_g)
236 31928 : ALLOCATE (phi_g)
237 31928 : CALL pw_pool%create_pw(phi_g)
238 : CALL pw_poisson_solve(poisson_env, rhob_g, vg_coulomb, phi_g, dphi_g, &
239 31928 : h_stress)
240 : !---------- END OF ELECTROSTATIC CALCULATION --------
241 :
242 : ! Atomic Energy
243 31928 : IF (atprop%energy) THEN
244 4872 : ALLOCATE (rpot)
245 232 : CALL rs_grid_create(rpot, rs_desc)
246 232 : CALL rs_grid_set_box(grid_spme, rs=rpot)
247 232 : CALL pw_multiply_with(phi_g, green%p3m_charge)
248 232 : CALL pw_transfer(phi_g, rhob_r)
249 232 : CALL transfer_pw2rs(rpot, rhob_r)
250 232 : ipart = 0
251 : DO
252 2335 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
253 2335 : IF (p1 == 0) EXIT
254 2103 : do_shell = (particle_set(p1)%shell_index /= 0)
255 2103 : IF (do_shell) CYCLE
256 : ! calculate function on small boxes
257 : CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
258 1023 : is_shell=.FALSE., unit_charge=.FALSE., charges=charges)
259 : ! integrate box and potential
260 1023 : CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fat1)
261 1255 : IF (atprop%energy) THEN
262 1023 : atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols
263 : END IF
264 : END DO
265 : ! Core-shell model
266 232 : IF (PRESENT(shell_particle_set) .AND. PRESENT(core_particle_set)) THEN
267 30 : ipart = 0
268 1110 : DO
269 1110 : CALL set_list(shell_particle_set, nshell, shell_center, p1_shell, rden, ipart)
270 1110 : IF (p1_shell == 0) EXIT
271 : CALL get_patch(shell_particle_set, shell_delta, green, p1_shell, rhos, &
272 1080 : is_core=.FALSE., is_shell=.TRUE., unit_charge=.FALSE.)
273 1080 : CALL dg_sum_patch_force_1d(rpot, rhos, shell_center(:, p1_shell), fat1)
274 1080 : p1 = shell_particle_set(p1_shell)%atom_index
275 1110 : IF (atprop%energy) THEN
276 1080 : atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols
277 : END IF
278 : END DO
279 30 : ipart = 0
280 1110 : DO
281 1110 : CALL set_list(core_particle_set, nshell, core_center, p1_shell, rden, ipart)
282 1110 : IF (p1_shell == 0) EXIT
283 : CALL get_patch(core_particle_set, core_delta, green, p1_shell, rhos, &
284 1080 : is_core=.TRUE., is_shell=.FALSE., unit_charge=.FALSE.)
285 1080 : CALL dg_sum_patch_force_1d(rpot, rhos, core_center(:, p1_shell), fat1)
286 1080 : p1 = core_particle_set(p1_shell)%atom_index
287 1110 : IF (atprop%energy) THEN
288 1080 : atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols
289 : END IF
290 : END DO
291 : END IF
292 232 : CALL rs_grid_release(rpot)
293 232 : DEALLOCATE (rpot)
294 : END IF
295 :
296 31928 : CALL pw_pool%give_back_pw(phi_g)
297 31928 : CALL pw_pool%give_back_pw(rhob_g)
298 31928 : DEALLOCATE (phi_g, rhob_g)
299 :
300 : !------------- STRESS TENSOR CALCULATION ------------
301 31928 : IF (use_virial) THEN
302 63560 : DO i = 1, 3
303 158900 : DO j = i, 3
304 95340 : f_stress(i, j) = pw_integral_a2b(dphi_g(i), dphi_g(j))
305 143010 : f_stress(j, i) = f_stress(i, j)
306 : END DO
307 : END DO
308 15890 : ffa = (1.0_dp/fourpi)*(0.5_dp/alpha)**2
309 206570 : f_stress = -ffa*f_stress
310 206570 : pv_g = h_stress + f_stress
311 : END IF
312 : !--------END OF STRESS TENSOR CALCULATION -----------
313 : ! move derivative of potential to real space grid and
314 : ! multiply by charge function in g-space
315 127712 : DO i = 1, 3
316 95784 : CALL rs_grid_create(drpot(i), rs_desc)
317 95784 : CALL rs_grid_set_box(grid_spme, rs=drpot(i))
318 95784 : CALL pw_multiply_with(dphi_g(i), green%p3m_charge)
319 95784 : CALL pw_transfer(dphi_g(i), rhob_r)
320 95784 : CALL pw_pool%give_back_pw(dphi_g(i))
321 127712 : CALL transfer_pw2rs(drpot(i), rhob_r)
322 : END DO
323 :
324 31928 : CALL pw_pool%give_back_pw(rhob_r)
325 31928 : DEALLOCATE (rhob_r)
326 : !----------------- FORCE CALCULATION ----------------
327 : ! initialize the forces
328 25832856 : fg_coulomb = 0.0_dp
329 : ! Particles
330 31928 : ipart = 0
331 : DO
332 3618054 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
333 3618054 : IF (p1 == 0) EXIT
334 :
335 3586126 : do_shell = (particle_set(p1)%shell_index /= 0)
336 3586126 : IF (do_shell) CYCLE
337 : ! calculate function on small boxes
338 : CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
339 3156071 : is_shell=.FALSE., unit_charge=.FALSE., charges=charges)
340 :
341 : ! add boxes to real space grid (big box)
342 3156071 : CALL dg_sum_patch_force_3d(drpot, rhos, center(:, p1), fat)
343 3156071 : fg_coulomb(1, p1) = fg_coulomb(1, p1) - fat(1)*dvols
344 3156071 : fg_coulomb(2, p1) = fg_coulomb(2, p1) - fat(2)*dvols
345 3586126 : fg_coulomb(3, p1) = fg_coulomb(3, p1) - fat(3)*dvols
346 : END DO
347 : ! Shell-Model
348 31928 : IF (PRESENT(shell_particle_set) .AND. (PRESENT(core_particle_set))) THEN
349 11732 : IF (PRESENT(fgshell_coulomb)) THEN
350 11732 : ipart = 0
351 3174012 : fgshell_coulomb = 0.0_dp
352 441787 : DO
353 : ! calculate function on small boxes
354 : CALL set_list(shell_particle_set, nshell, shell_center, p1_shell, &
355 441787 : rden, ipart)
356 441787 : IF (p1_shell == 0) EXIT
357 :
358 : CALL get_patch(shell_particle_set, shell_delta, green, &
359 430055 : p1_shell, rhos, is_core=.FALSE., is_shell=.TRUE., unit_charge=.FALSE.)
360 :
361 : ! add boxes to real space grid (big box)
362 430055 : CALL dg_sum_patch_force_3d(drpot, rhos, shell_center(:, p1_shell), fat)
363 430055 : fgshell_coulomb(1, p1_shell) = fgshell_coulomb(1, p1_shell) - fat(1)*dvols
364 430055 : fgshell_coulomb(2, p1_shell) = fgshell_coulomb(2, p1_shell) - fat(2)*dvols
365 430055 : fgshell_coulomb(3, p1_shell) = fgshell_coulomb(3, p1_shell) - fat(3)*dvols
366 :
367 : END DO
368 : END IF
369 11732 : IF (PRESENT(fgcore_coulomb)) THEN
370 11732 : ipart = 0
371 3174012 : fgcore_coulomb = 0.0_dp
372 441787 : DO
373 : ! calculate function on small boxes
374 : CALL set_list(core_particle_set, nshell, core_center, p1_shell, &
375 441787 : rden, ipart)
376 441787 : IF (p1_shell == 0) EXIT
377 :
378 : CALL get_patch(core_particle_set, core_delta, green, &
379 430055 : p1_shell, rhos, is_core=.TRUE., is_shell=.FALSE., unit_charge=.FALSE.)
380 :
381 : ! add boxes to real space grid (big box)
382 430055 : CALL dg_sum_patch_force_3d(drpot, rhos, core_center(:, p1_shell), fat)
383 430055 : fgcore_coulomb(1, p1_shell) = fgcore_coulomb(1, p1_shell) - fat(1)*dvols
384 430055 : fgcore_coulomb(2, p1_shell) = fgcore_coulomb(2, p1_shell) - fat(2)*dvols
385 430055 : fgcore_coulomb(3, p1_shell) = fgcore_coulomb(3, p1_shell) - fat(3)*dvols
386 : END DO
387 : END IF
388 :
389 : END IF
390 : !--------------END OF FORCE CALCULATION -------------
391 :
392 : !------------------CLEANING UP ----------------------
393 31928 : CALL rs_grid_release(rden)
394 31928 : DEALLOCATE (rden)
395 127712 : DO i = 1, 3
396 127712 : CALL rs_grid_release(drpot(i))
397 : END DO
398 :
399 31928 : DEALLOCATE (rhos)
400 31928 : DEALLOCATE (center, delta)
401 31928 : IF (ALLOCATED(shell_center)) THEN
402 11732 : DEALLOCATE (shell_center, shell_delta)
403 : END IF
404 31928 : IF (ALLOCATED(core_center)) THEN
405 11732 : DEALLOCATE (core_center, core_delta)
406 : END IF
407 31928 : CALL timestop(handle)
408 :
409 159640 : END SUBROUTINE spme_evaluate
410 :
411 : ! **************************************************************************************************
412 : !> \brief Calculate the electrostatic potential from particles A (charge A) at
413 : !> positions of particles B
414 : !> \param ewald_env ...
415 : !> \param ewald_pw ...
416 : !> \param box ...
417 : !> \param particle_set_a ...
418 : !> \param charges_a ...
419 : !> \param particle_set_b ...
420 : !> \param potential ...
421 : !> \par History
422 : !> JGH (23-Oct-2014) : adapted from SPME evaluate
423 : !> \author JGH (23-Oct-2014)
424 : ! **************************************************************************************************
425 13942 : SUBROUTINE spme_potential(ewald_env, ewald_pw, box, particle_set_a, charges_a, &
426 13942 : particle_set_b, potential)
427 :
428 : TYPE(ewald_environment_type), POINTER :: ewald_env
429 : TYPE(ewald_pw_type), POINTER :: ewald_pw
430 : TYPE(cell_type), POINTER :: box
431 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set_a
432 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), POINTER :: charges_a
433 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set_b
434 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: potential
435 :
436 : CHARACTER(len=*), PARAMETER :: routineN = 'spme_potential'
437 :
438 : INTEGER :: handle, ipart, n, npart_a, npart_b, &
439 : o_spline, p1
440 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: center
441 : INTEGER, DIMENSION(3) :: npts
442 : REAL(KIND=dp) :: alpha, dvols, fat1
443 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: delta
444 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rhos
445 : TYPE(greens_fn_type), POINTER :: green
446 : TYPE(mp_comm_type) :: group
447 : TYPE(pw_c1d_gs_type), POINTER :: phi_g, rhob_g
448 : TYPE(pw_grid_type), POINTER :: grid_spme
449 : TYPE(pw_poisson_type), POINTER :: poisson_env
450 : TYPE(pw_pool_type), POINTER :: pw_pool
451 : TYPE(pw_r3d_rs_type), POINTER :: rhob_r
452 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
453 : TYPE(realspace_grid_type), POINTER :: rden, rpot
454 :
455 13942 : NULLIFY (grid_spme, green, poisson_env, phi_g, rhob_g, rhob_r, pw_pool, &
456 13942 : rden, rpot)
457 13942 : CALL timeset(routineN, handle)
458 13942 : CALL cite_reference(Essmann1995)
459 :
460 : !-------------- INITIALISATION ---------------------
461 13942 : CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group)
462 13942 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, poisson_env=poisson_env)
463 13942 : CALL pw_poisson_rebuild(poisson_env)
464 13942 : green => poisson_env%green_fft
465 13942 : grid_spme => pw_pool%pw_grid
466 :
467 13942 : npart_a = SIZE(particle_set_a)
468 13942 : npart_b = SIZE(particle_set_b)
469 :
470 13942 : CALL get_pw_grid_info(grid_spme, npts=npts, dvol=dvols)
471 :
472 13942 : n = o_spline
473 69710 : ALLOCATE (rhos(n, n, n))
474 292782 : ALLOCATE (rden)
475 13942 : CALL rs_grid_create(rden, rs_desc)
476 13942 : CALL rs_grid_set_box(grid_spme, rs=rden)
477 13942 : CALL rs_grid_zero(rden)
478 :
479 69710 : ALLOCATE (center(3, npart_a), delta(3, npart_a))
480 13942 : CALL get_center(particle_set_a, box, center, delta, npts, n)
481 :
482 : !-------------- DENSITY CALCULATION ----------------
483 13942 : ipart = 0
484 : ! Particles
485 43272 : DO
486 57214 : CALL set_list(particle_set_a, npart_a, center, p1, rden, ipart)
487 57214 : IF (p1 == 0) EXIT
488 :
489 : ! calculate function on small boxes
490 43272 : CALL get_patch(particle_set_a, delta, green, p1, rhos, charges_a)
491 :
492 : ! add boxes to real space grid (big box)
493 57214 : CALL dg_sum_patch(rden, rhos, center(:, p1))
494 : END DO
495 13942 : DEALLOCATE (center, delta)
496 : !----------- END OF DENSITY CALCULATION -------------
497 :
498 13942 : NULLIFY (rhob_r)
499 13942 : ALLOCATE (rhob_r)
500 13942 : CALL pw_pool%create_pw(rhob_r)
501 13942 : CALL transfer_rs2pw(rden, rhob_r)
502 : ! transform density to G space and add charge function
503 13942 : NULLIFY (rhob_g)
504 13942 : ALLOCATE (rhob_g)
505 13942 : CALL pw_pool%create_pw(rhob_g)
506 13942 : CALL pw_transfer(rhob_r, rhob_g)
507 : ! update charge function
508 13942 : CALL pw_multiply_with(rhob_g, green%p3m_charge)
509 :
510 : !-------------- ELECTROSTATIC CALCULATION -----------
511 : ! allocate intermediate arrays
512 13942 : NULLIFY (phi_g)
513 13942 : ALLOCATE (phi_g)
514 13942 : CALL pw_pool%create_pw(phi_g)
515 13942 : CALL pw_poisson_solve(poisson_env, density=rhob_g, vhartree=phi_g)
516 : !---------- END OF ELECTROSTATIC CALCULATION --------
517 :
518 : !-------------- POTENTAIL AT ATOMIC POSITIONS -------
519 69710 : ALLOCATE (center(3, npart_b), delta(3, npart_b))
520 13942 : CALL get_center(particle_set_b, box, center, delta, npts, n)
521 13942 : rpot => rden
522 13942 : CALL pw_multiply_with(phi_g, green%p3m_charge)
523 13942 : CALL pw_transfer(phi_g, rhob_r)
524 13942 : CALL transfer_pw2rs(rpot, rhob_r)
525 13942 : ipart = 0
526 43272 : DO
527 57214 : CALL set_list(particle_set_b, npart_b, center, p1, rden, ipart)
528 57214 : IF (p1 == 0) EXIT
529 : ! calculate function on small boxes
530 : CALL get_patch(particle_set_b, delta, green, p1, rhos, &
531 43272 : is_core=.FALSE., is_shell=.FALSE., unit_charge=.TRUE.)
532 : ! integrate box and potential
533 43272 : CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fat1)
534 43272 : potential(p1) = potential(p1) + fat1*dvols
535 : END DO
536 :
537 : !------------------CLEANING UP ----------------------
538 13942 : CALL pw_pool%give_back_pw(phi_g)
539 13942 : CALL pw_pool%give_back_pw(rhob_g)
540 13942 : CALL pw_pool%give_back_pw(rhob_r)
541 13942 : DEALLOCATE (phi_g, rhob_g, rhob_r)
542 13942 : CALL rs_grid_release(rden)
543 13942 : DEALLOCATE (rden)
544 :
545 13942 : DEALLOCATE (rhos)
546 13942 : DEALLOCATE (center, delta)
547 13942 : CALL timestop(handle)
548 :
549 41826 : END SUBROUTINE spme_potential
550 :
551 : ! **************************************************************************************************
552 : !> \brief Calculate the forces on particles B for the electrostatic interaction
553 : !> betrween particles A and B
554 : !> \param ewald_env ...
555 : !> \param ewald_pw ...
556 : !> \param box ...
557 : !> \param particle_set_a ...
558 : !> \param charges_a ...
559 : !> \param particle_set_b ...
560 : !> \param charges_b ...
561 : !> \param forces_b ...
562 : !> \par History
563 : !> JGH (27-Oct-2014) : adapted from SPME evaluate
564 : !> \author JGH (23-Oct-2014)
565 : ! **************************************************************************************************
566 160 : SUBROUTINE spme_forces(ewald_env, ewald_pw, box, particle_set_a, charges_a, &
567 160 : particle_set_b, charges_b, forces_b)
568 :
569 : TYPE(ewald_environment_type), POINTER :: ewald_env
570 : TYPE(ewald_pw_type), POINTER :: ewald_pw
571 : TYPE(cell_type), POINTER :: box
572 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set_a
573 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), POINTER :: charges_a
574 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set_b
575 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), POINTER :: charges_b
576 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: forces_b
577 :
578 : CHARACTER(len=*), PARAMETER :: routineN = 'spme_forces'
579 :
580 : INTEGER :: handle, i, ipart, n, npart_a, npart_b, &
581 : o_spline, p1
582 160 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: center
583 : INTEGER, DIMENSION(3) :: npts
584 : REAL(KIND=dp) :: alpha, dvols
585 160 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: delta
586 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rhos
587 : REAL(KIND=dp), DIMENSION(3) :: fat
588 : TYPE(greens_fn_type), POINTER :: green
589 : TYPE(mp_comm_type) :: group
590 640 : TYPE(pw_c1d_gs_type), DIMENSION(3) :: dphi_g
591 : TYPE(pw_c1d_gs_type), POINTER :: phi_g, rhob_g
592 : TYPE(pw_grid_type), POINTER :: grid_spme
593 : TYPE(pw_poisson_type), POINTER :: poisson_env
594 : TYPE(pw_pool_type), POINTER :: pw_pool
595 : TYPE(pw_r3d_rs_type), POINTER :: rhob_r
596 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
597 5760 : TYPE(realspace_grid_type), DIMENSION(3) :: drpot
598 : TYPE(realspace_grid_type), POINTER :: rden
599 :
600 160 : NULLIFY (grid_spme, green, poisson_env, phi_g, rhob_g, rhob_r, &
601 160 : pw_pool, rden)
602 160 : CALL timeset(routineN, handle)
603 160 : CALL cite_reference(Essmann1995)
604 :
605 : !-------------- INITIALISATION ---------------------
606 160 : CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group)
607 160 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, poisson_env=poisson_env)
608 160 : CALL pw_poisson_rebuild(poisson_env)
609 160 : green => poisson_env%green_fft
610 160 : grid_spme => pw_pool%pw_grid
611 :
612 160 : npart_a = SIZE(particle_set_a)
613 160 : npart_b = SIZE(particle_set_b)
614 :
615 160 : CALL get_pw_grid_info(grid_spme, npts=npts, dvol=dvols)
616 :
617 160 : n = o_spline
618 800 : ALLOCATE (rhos(n, n, n))
619 3360 : ALLOCATE (rden)
620 160 : CALL rs_grid_create(rden, rs_desc)
621 160 : CALL rs_grid_set_box(grid_spme, rs=rden)
622 160 : CALL rs_grid_zero(rden)
623 :
624 800 : ALLOCATE (center(3, npart_a), delta(3, npart_a))
625 160 : CALL get_center(particle_set_a, box, center, delta, npts, n)
626 :
627 : !-------------- DENSITY CALCULATION ----------------
628 160 : ipart = 0
629 : ! Particles
630 380 : DO
631 540 : CALL set_list(particle_set_a, npart_a, center, p1, rden, ipart)
632 540 : IF (p1 == 0) EXIT
633 :
634 : ! calculate function on small boxes
635 380 : CALL get_patch(particle_set_a, delta, green, p1, rhos, charges_a)
636 :
637 : ! add boxes to real space grid (big box)
638 540 : CALL dg_sum_patch(rden, rhos, center(:, p1))
639 : END DO
640 160 : DEALLOCATE (center, delta)
641 : !----------- END OF DENSITY CALCULATION -------------
642 :
643 160 : NULLIFY (rhob_r)
644 160 : ALLOCATE (rhob_r)
645 160 : CALL pw_pool%create_pw(rhob_r)
646 160 : CALL transfer_rs2pw(rden, rhob_r)
647 : ! transform density to G space and add charge function
648 160 : NULLIFY (rhob_g)
649 160 : ALLOCATE (rhob_g)
650 160 : CALL pw_pool%create_pw(rhob_g)
651 160 : CALL pw_transfer(rhob_r, rhob_g)
652 : ! update charge function
653 160 : CALL pw_multiply_with(rhob_g, green%p3m_charge)
654 :
655 : !-------------- ELECTROSTATIC CALCULATION -----------
656 : ! allocate intermediate arrays
657 640 : DO i = 1, 3
658 640 : CALL pw_pool%create_pw(dphi_g(i))
659 : END DO
660 160 : NULLIFY (phi_g)
661 160 : ALLOCATE (phi_g)
662 160 : CALL pw_pool%create_pw(phi_g)
663 : CALL pw_poisson_solve(poisson_env, density=rhob_g, vhartree=phi_g, &
664 160 : dvhartree=dphi_g)
665 : !---------- END OF ELECTROSTATIC CALCULATION --------
666 : ! move derivative of potential to real space grid and
667 : ! multiply by charge function in g-space
668 640 : DO i = 1, 3
669 480 : CALL rs_grid_create(drpot(i), rs_desc)
670 480 : CALL rs_grid_set_box(grid_spme, rs=drpot(i))
671 480 : CALL pw_multiply_with(dphi_g(i), green%p3m_charge)
672 480 : CALL pw_transfer(dphi_g(i), rhob_r)
673 480 : CALL pw_pool%give_back_pw(dphi_g(i))
674 640 : CALL transfer_pw2rs(drpot(i), rhob_r)
675 : END DO
676 : !----------------- FORCE CALCULATION ----------------
677 800 : ALLOCATE (center(3, npart_b), delta(3, npart_b))
678 160 : CALL get_center(particle_set_b, box, center, delta, npts, n)
679 160 : ipart = 0
680 380 : DO
681 540 : CALL set_list(particle_set_b, npart_b, center, p1, rden, ipart)
682 540 : IF (p1 == 0) EXIT
683 : ! calculate function on small boxes
684 : CALL get_patch(particle_set_b, delta, green, p1, rhos, &
685 380 : is_core=.FALSE., is_shell=.FALSE., unit_charge=.FALSE., charges=charges_b)
686 : ! add boxes to real space grid (big box)
687 380 : CALL dg_sum_patch_force_3d(drpot, rhos, center(:, p1), fat)
688 380 : forces_b(1, p1) = forces_b(1, p1) - fat(1)*dvols
689 380 : forces_b(2, p1) = forces_b(2, p1) - fat(2)*dvols
690 540 : forces_b(3, p1) = forces_b(3, p1) - fat(3)*dvols
691 : END DO
692 : !------------------CLEANING UP ----------------------
693 640 : DO i = 1, 3
694 640 : CALL rs_grid_release(drpot(i))
695 : END DO
696 160 : CALL pw_pool%give_back_pw(phi_g)
697 160 : CALL pw_pool%give_back_pw(rhob_g)
698 160 : CALL pw_pool%give_back_pw(rhob_r)
699 160 : DEALLOCATE (phi_g, rhob_g, rhob_r)
700 160 : CALL rs_grid_release(rden)
701 160 : DEALLOCATE (rden)
702 :
703 160 : DEALLOCATE (rhos)
704 160 : DEALLOCATE (center, delta)
705 160 : CALL timestop(handle)
706 :
707 2400 : END SUBROUTINE spme_forces
708 :
709 : ! **************************************************************************************************
710 : !> \brief Internal Virial for 1/2 [rho||rho] (rho=mcharge)
711 : !> \param ewald_env ...
712 : !> \param ewald_pw ...
713 : !> \param particle_set ...
714 : !> \param box ...
715 : !> \param mcharge ...
716 : !> \param virial ...
717 : ! **************************************************************************************************
718 52 : SUBROUTINE spme_virial(ewald_env, ewald_pw, particle_set, box, mcharge, virial)
719 :
720 : TYPE(ewald_environment_type), POINTER :: ewald_env
721 : TYPE(ewald_pw_type), POINTER :: ewald_pw
722 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
723 : TYPE(cell_type), POINTER :: box
724 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: mcharge
725 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: virial
726 :
727 : CHARACTER(len=*), PARAMETER :: routineN = 'spme_virial'
728 :
729 : INTEGER :: handle, i, ipart, j, n, npart, o_spline, &
730 : p1
731 52 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: center
732 : INTEGER, DIMENSION(3) :: npts
733 : REAL(KIND=dp) :: alpha, dvols, ffa, vgc
734 52 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: delta
735 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rhos
736 : REAL(KIND=dp), DIMENSION(3, 3) :: f_stress, h_stress
737 : TYPE(greens_fn_type), POINTER :: green
738 : TYPE(mp_comm_type) :: group
739 : TYPE(mp_para_env_type), POINTER :: para_env
740 208 : TYPE(pw_c1d_gs_type), DIMENSION(3) :: dphi_g
741 : TYPE(pw_c1d_gs_type), POINTER :: phi_g, rhob_g
742 : TYPE(pw_grid_type), POINTER :: grid_spme
743 : TYPE(pw_poisson_type), POINTER :: poisson_env
744 : TYPE(pw_pool_type), POINTER :: pw_pool
745 : TYPE(pw_r3d_rs_type), POINTER :: rhob_r
746 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
747 2028 : TYPE(realspace_grid_type) :: rden, rpot
748 :
749 52 : CALL timeset(routineN, handle)
750 : !-------------- INITIALISATION ---------------------
751 52 : virial = 0.0_dp
752 : CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group, &
753 52 : para_env=para_env)
754 52 : NULLIFY (green, poisson_env, pw_pool)
755 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, &
756 52 : poisson_env=poisson_env)
757 52 : CALL pw_poisson_rebuild(poisson_env)
758 52 : green => poisson_env%green_fft
759 52 : grid_spme => pw_pool%pw_grid
760 :
761 52 : CALL get_pw_grid_info(grid_spme, dvol=dvols, npts=npts)
762 :
763 52 : npart = SIZE(particle_set)
764 :
765 52 : n = o_spline
766 260 : ALLOCATE (rhos(n, n, n))
767 :
768 52 : CALL rs_grid_create(rden, rs_desc)
769 52 : CALL rs_grid_set_box(grid_spme, rs=rden)
770 52 : CALL rs_grid_zero(rden)
771 :
772 260 : ALLOCATE (center(3, npart), delta(3, npart))
773 52 : CALL get_center(particle_set, box, center, delta, npts, n)
774 :
775 : !-------------- DENSITY CALCULATION ----------------
776 52 : ipart = 0
777 194 : DO
778 246 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
779 246 : IF (p1 == 0) EXIT
780 :
781 : ! calculate function on small boxes
782 : CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
783 194 : is_shell=.FALSE., unit_charge=.TRUE.)
784 16490 : rhos(:, :, :) = rhos(:, :, :)*mcharge(p1)
785 :
786 : ! add boxes to real space grid (big box)
787 246 : CALL dg_sum_patch(rden, rhos, center(:, p1))
788 : END DO
789 :
790 52 : NULLIFY (rhob_r)
791 52 : ALLOCATE (rhob_r)
792 52 : CALL pw_pool%create_pw(rhob_r)
793 :
794 52 : CALL transfer_rs2pw(rden, rhob_r)
795 :
796 : ! transform density to G space and add charge function
797 52 : NULLIFY (rhob_g)
798 52 : ALLOCATE (rhob_g)
799 52 : CALL pw_pool%create_pw(rhob_g)
800 52 : CALL pw_transfer(rhob_r, rhob_g)
801 : ! update charge function
802 52 : CALL pw_multiply_with(rhob_g, green%p3m_charge)
803 :
804 : !-------------- ELECTROSTATIC CALCULATION -----------
805 :
806 : ! allocate intermediate arrays
807 208 : DO i = 1, 3
808 208 : CALL pw_pool%create_pw(dphi_g(i))
809 : END DO
810 52 : NULLIFY (phi_g)
811 52 : ALLOCATE (phi_g)
812 52 : CALL pw_pool%create_pw(phi_g)
813 52 : CALL pw_poisson_solve(poisson_env, rhob_g, vgc, phi_g, dphi_g, h_stress=h_stress)
814 :
815 52 : CALL rs_grid_create(rpot, rs_desc)
816 52 : CALL rs_grid_set_box(grid_spme, rs=rpot)
817 :
818 52 : CALL pw_pool%give_back_pw(rhob_g)
819 52 : DEALLOCATE (rhob_g)
820 :
821 52 : CALL rs_grid_zero(rpot)
822 52 : CALL pw_multiply_with(phi_g, green%p3m_charge)
823 52 : CALL pw_transfer(phi_g, rhob_r)
824 52 : CALL pw_pool%give_back_pw(phi_g)
825 52 : DEALLOCATE (phi_g)
826 52 : CALL transfer_pw2rs(rpot, rhob_r)
827 :
828 : !---------- END OF ELECTROSTATIC CALCULATION --------
829 :
830 : !------------- STRESS TENSOR CALCULATION ------------
831 :
832 208 : DO i = 1, 3
833 520 : DO j = i, 3
834 312 : f_stress(i, j) = pw_integral_a2b(dphi_g(i), dphi_g(j))
835 468 : f_stress(j, i) = f_stress(i, j)
836 : END DO
837 : END DO
838 52 : ffa = (1.0_dp/fourpi)*(0.5_dp/alpha)**2
839 676 : virial = virial - (ffa*f_stress - h_stress)/REAL(para_env%num_pe, dp)
840 :
841 : !--------END OF STRESS TENSOR CALCULATION -----------
842 :
843 208 : DO i = 1, 3
844 208 : CALL pw_pool%give_back_pw(dphi_g(i))
845 : END DO
846 52 : CALL pw_pool%give_back_pw(rhob_r)
847 52 : DEALLOCATE (rhob_r)
848 :
849 52 : CALL rs_grid_release(rden)
850 52 : CALL rs_grid_release(rpot)
851 52 : DEALLOCATE (rhos)
852 52 : DEALLOCATE (center, delta)
853 :
854 52 : CALL timestop(handle)
855 :
856 156 : END SUBROUTINE spme_virial
857 :
858 : ! **************************************************************************************************
859 : !> \brief Calculates local density in a small box
860 : !> \param part ...
861 : !> \param delta ...
862 : !> \param green ...
863 : !> \param p ...
864 : !> \param rhos ...
865 : !> \param is_core ...
866 : !> \param is_shell ...
867 : !> \param unit_charge ...
868 : !> \param charges ...
869 : !> \par History
870 : !> none
871 : !> \author JGH (21-Mar-2001)
872 : ! **************************************************************************************************
873 8371189 : SUBROUTINE get_patch_a(part, delta, green, p, rhos, is_core, is_shell, &
874 : unit_charge, charges)
875 :
876 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: part
877 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: delta
878 : TYPE(greens_fn_type), INTENT(IN) :: green
879 : INTEGER, INTENT(IN) :: p
880 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: rhos
881 : LOGICAL, INTENT(IN) :: is_core, is_shell, unit_charge
882 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges
883 :
884 : INTEGER :: nbox
885 : LOGICAL :: use_charge_array
886 : REAL(KIND=dp) :: q
887 : REAL(KIND=dp), DIMENSION(3) :: r
888 : TYPE(shell_kind_type), POINTER :: shell
889 :
890 8371189 : NULLIFY (shell)
891 8371189 : use_charge_array = .FALSE.
892 8371189 : IF (PRESENT(charges)) use_charge_array = ASSOCIATED(charges)
893 8371189 : IF (is_core .AND. is_shell) THEN
894 0 : CPABORT("Shell-model: cannot be core and shell simultaneously")
895 : END IF
896 :
897 8371189 : nbox = SIZE(rhos, 1)
898 33484756 : r = part(p)%r
899 8371189 : q = 1.0_dp
900 8371189 : IF (.NOT. unit_charge) THEN
901 8035925 : IF (is_core) THEN
902 861190 : CALL get_atomic_kind(atomic_kind=part(p)%atomic_kind, shell=shell)
903 861190 : q = shell%charge_core
904 7174735 : ELSE IF (is_shell) THEN
905 861190 : CALL get_atomic_kind(atomic_kind=part(p)%atomic_kind, shell=shell)
906 861190 : q = shell%charge_shell
907 : ELSE
908 6313545 : CALL get_atomic_kind(atomic_kind=part(p)%atomic_kind, qeff=q)
909 : END IF
910 8035925 : IF (use_charge_array) q = charges(p)
911 : END IF
912 8371189 : CALL spme_get_patch(rhos, nbox, delta(:, p), q, green%p3m_coeff)
913 :
914 8371189 : END SUBROUTINE get_patch_a
915 :
916 : ! **************************************************************************************************
917 : !> \brief Calculates local density in a small box
918 : !> \param part ...
919 : !> \param delta ...
920 : !> \param green ...
921 : !> \param p ...
922 : !> \param rhos ...
923 : !> \param charges ...
924 : !> \par History
925 : !> none
926 : !> \author JGH (21-Mar-2001)
927 : ! **************************************************************************************************
928 43652 : SUBROUTINE get_patch_b(part, delta, green, p, rhos, charges)
929 :
930 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: part
931 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: delta
932 : TYPE(greens_fn_type), INTENT(IN) :: green
933 : INTEGER, INTENT(IN) :: p
934 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: rhos
935 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges
936 :
937 : INTEGER :: nbox
938 : REAL(KIND=dp) :: q
939 : REAL(KIND=dp), DIMENSION(3) :: r
940 :
941 43652 : CPASSERT(ASSOCIATED(charges))
942 43652 : nbox = SIZE(rhos, 1)
943 174608 : r = part(p)%r
944 43652 : q = charges(p)
945 43652 : CALL spme_get_patch(rhos, nbox, delta(:, p), q, green%p3m_coeff)
946 :
947 43652 : END SUBROUTINE get_patch_b
948 :
949 : ! **************************************************************************************************
950 : !> \brief Calculates SPME charge assignment
951 : !> \param rhos ...
952 : !> \param n ...
953 : !> \param delta ...
954 : !> \param q ...
955 : !> \param coeff ...
956 : !> \par History
957 : !> DG (29-Mar-2001) : code implemented
958 : !> \author JGH (22-Mar-2001)
959 : ! **************************************************************************************************
960 8414841 : SUBROUTINE spme_get_patch(rhos, n, delta, q, coeff)
961 :
962 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: rhos
963 : INTEGER, INTENT(IN) :: n
964 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: delta
965 : REAL(KIND=dp), INTENT(IN) :: q
966 : REAL(KIND=dp), DIMENSION(-(n-1):n-1, 0:n-1), &
967 : INTENT(IN) :: coeff
968 :
969 : INTEGER, PARAMETER :: nmax = 12
970 :
971 : INTEGER :: i, i1, i2, i3, j, l
972 : REAL(KIND=dp) :: r2, r3
973 : REAL(KIND=dp), DIMENSION(3, -nmax:nmax) :: w_assign
974 : REAL(KIND=dp), DIMENSION(3, 0:nmax-1) :: deltal
975 : REAL(KIND=dp), DIMENSION(3, 1:nmax) :: f_assign
976 :
977 8414841 : IF (n > nmax) THEN
978 0 : CPABORT("nmax value too small")
979 : END IF
980 : ! calculate the assignment function values and
981 : ! the charges on the grid (small box)
982 :
983 8414841 : deltal(1, 0) = 1.0_dp
984 8414841 : deltal(2, 0) = 1.0_dp
985 8414841 : deltal(3, 0) = 1.0_dp
986 44362920 : DO l = 1, n - 1
987 35948079 : deltal(1, l) = deltal(1, l - 1)*delta(1)
988 35948079 : deltal(2, l) = deltal(2, l - 1)*delta(2)
989 44362920 : deltal(3, l) = deltal(3, l - 1)*delta(3)
990 : END DO
991 :
992 8414841 : w_assign = 0.0_dp
993 52777761 : DO j = -(n - 1), n - 1, 2
994 294542195 : DO l = 0, n - 1
995 241764434 : w_assign(1, j) = w_assign(1, j) + coeff(j, l)*deltal(1, l)
996 241764434 : w_assign(2, j) = w_assign(2, j) + coeff(j, l)*deltal(2, l)
997 286127354 : w_assign(3, j) = w_assign(3, j) + coeff(j, l)*deltal(3, l)
998 : END DO
999 : END DO
1000 52777761 : DO i = 1, n
1001 44362920 : j = n + 1 - 2*i
1002 44362920 : f_assign(1, i) = w_assign(1, j)
1003 44362920 : f_assign(2, i) = w_assign(2, j)
1004 52777761 : f_assign(3, i) = w_assign(3, j)
1005 : END DO
1006 :
1007 52777761 : DO i3 = 1, n
1008 44362920 : r3 = q*f_assign(3, i3)
1009 294542195 : DO i2 = 1, n
1010 241764434 : r2 = r3*f_assign(2, i2)
1011 1640002544 : DO i1 = 1, n
1012 1595639624 : rhos(i1, i2, i3) = r2*f_assign(1, i1)
1013 : END DO
1014 : END DO
1015 : END DO
1016 :
1017 8414841 : END SUBROUTINE spme_get_patch
1018 :
1019 : END MODULE spme
1020 :
|
