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 28842 : SUBROUTINE spme_evaluate(ewald_env, ewald_pw, box, particle_set, &
92 28842 : fg_coulomb, vg_coulomb, pv_g, shell_particle_set, core_particle_set, &
93 28842 : 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 57684 : 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 28842 : 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 115368 : 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 1038312 : TYPE(realspace_grid_type), DIMENSION(3) :: drpot
132 : TYPE(realspace_grid_type), POINTER :: rden, rpot
133 :
134 28842 : NULLIFY (grid_spme, green, poisson_env, phi_g, rhob_g, rhob_r, pw_pool, rden, rpot)
135 28842 : CALL timeset(routineN, handle)
136 28842 : CALL cite_reference(Essmann1995)
137 :
138 : !-------------- INITIALISATION ---------------------
139 : CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, &
140 28842 : group=group)
141 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, &
142 28842 : poisson_env=poisson_env)
143 28842 : CALL pw_poisson_rebuild(poisson_env)
144 28842 : green => poisson_env%green_fft
145 28842 : grid_spme => pw_pool%pw_grid
146 :
147 28842 : npart = SIZE(particle_set)
148 :
149 28842 : CALL get_pw_grid_info(grid_spme, npts=npts, dvol=dvols)
150 :
151 28842 : n = o_spline
152 144210 : ALLOCATE (rhos(n, n, n))
153 605682 : ALLOCATE (rden)
154 28842 : CALL rs_grid_create(rden, rs_desc)
155 28842 : CALL rs_grid_set_box(grid_spme, rs=rden)
156 28842 : CALL rs_grid_zero(rden)
157 :
158 144210 : ALLOCATE (center(3, npart), delta(3, npart))
159 28842 : CALL get_center(particle_set, box, center, delta, npts, n)
160 :
161 28842 : IF (PRESENT(shell_particle_set) .AND. (PRESENT(core_particle_set))) THEN
162 10396 : CPASSERT(ASSOCIATED(shell_particle_set))
163 10396 : CPASSERT(ASSOCIATED(core_particle_set))
164 10396 : nshell = SIZE(shell_particle_set)
165 10396 : ncore = SIZE(core_particle_set)
166 10396 : CPASSERT(nshell == ncore)
167 51980 : ALLOCATE (shell_center(3, nshell), shell_delta(3, nshell))
168 10396 : CALL get_center(shell_particle_set, box, shell_center, shell_delta, npts, n)
169 31188 : ALLOCATE (core_center(3, nshell), core_delta(3, nshell))
170 10396 : CALL get_center(core_particle_set, box, core_center, core_delta, npts, n)
171 : END IF
172 :
173 : !-------------- DENSITY CALCULATION ----------------
174 28842 : ipart = 0
175 : ! Particles
176 : DO
177 3582406 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
178 3582406 : IF (p1 == 0) EXIT
179 :
180 3553564 : do_shell = (particle_set(p1)%shell_index /= 0)
181 3553564 : IF (do_shell) CYCLE
182 : ! calculate function on small boxes
183 : CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
184 3134675 : is_shell=.FALSE., unit_charge=.FALSE., charges=charges)
185 :
186 : ! add boxes to real space grid (big box)
187 3582406 : CALL dg_sum_patch(rden, rhos, center(:, p1))
188 : END DO
189 : ! Shell-Model
190 28842 : IF (PRESENT(shell_particle_set) .AND. PRESENT(core_particle_set)) THEN
191 10396 : ipart = 0
192 429285 : DO
193 : ! calculate function on small boxes
194 : CALL set_list(shell_particle_set, nshell, shell_center, p1_shell, &
195 429285 : rden, ipart)
196 429285 : IF (p1_shell == 0) EXIT
197 : CALL get_patch(shell_particle_set, shell_delta, green, p1_shell, rhos, &
198 418889 : is_core=.FALSE., is_shell=.TRUE., unit_charge=.FALSE.)
199 :
200 : ! add boxes to real space grid (big box)
201 429285 : CALL dg_sum_patch(rden, rhos, shell_center(:, p1_shell))
202 : END DO
203 10396 : ipart = 0
204 429285 : DO
205 : ! calculate function on small boxes
206 : CALL set_list(core_particle_set, nshell, core_center, p1_shell, &
207 429285 : rden, ipart)
208 429285 : IF (p1_shell == 0) EXIT
209 : CALL get_patch(core_particle_set, core_delta, green, p1_shell, rhos, &
210 418889 : is_core=.TRUE., is_shell=.FALSE., unit_charge=.FALSE.)
211 :
212 : ! add boxes to real space grid (big box)
213 429285 : CALL dg_sum_patch(rden, rhos, core_center(:, p1_shell))
214 : END DO
215 : END IF
216 : !----------- END OF DENSITY CALCULATION -------------
217 :
218 28842 : NULLIFY (rhob_r)
219 28842 : ALLOCATE (rhob_r)
220 28842 : CALL pw_pool%create_pw(rhob_r)
221 28842 : CALL transfer_rs2pw(rden, rhob_r)
222 : ! transform density to G space and add charge function
223 28842 : NULLIFY (rhob_g)
224 28842 : ALLOCATE (rhob_g)
225 28842 : CALL pw_pool%create_pw(rhob_g)
226 28842 : CALL pw_transfer(rhob_r, rhob_g)
227 : ! update charge function
228 28842 : CALL pw_multiply_with(rhob_g, green%p3m_charge)
229 :
230 : !-------------- ELECTROSTATIC CALCULATION -----------
231 : ! allocate intermediate arrays
232 115368 : DO i = 1, 3
233 115368 : CALL pw_pool%create_pw(dphi_g(i))
234 : END DO
235 28842 : NULLIFY (phi_g)
236 28842 : ALLOCATE (phi_g)
237 28842 : CALL pw_pool%create_pw(phi_g)
238 : CALL pw_poisson_solve(poisson_env, rhob_g, vg_coulomb, phi_g, dphi_g, &
239 28842 : h_stress)
240 : !---------- END OF ELECTROSTATIC CALCULATION --------
241 :
242 : ! Atomic Energy
243 28842 : 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 28842 : CALL pw_pool%give_back_pw(phi_g)
297 28842 : CALL pw_pool%give_back_pw(rhob_g)
298 28842 : DEALLOCATE (phi_g, rhob_g)
299 :
300 : !------------- STRESS TENSOR CALCULATION ------------
301 28842 : IF (use_virial) THEN
302 49192 : DO i = 1, 3
303 122980 : DO j = i, 3
304 73788 : f_stress(i, j) = pw_integral_a2b(dphi_g(i), dphi_g(j))
305 110682 : f_stress(j, i) = f_stress(i, j)
306 : END DO
307 : END DO
308 12298 : ffa = (1.0_dp/fourpi)*(0.5_dp/alpha)**2
309 159874 : f_stress = -ffa*f_stress
310 159874 : 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 115368 : DO i = 1, 3
316 86526 : CALL rs_grid_create(drpot(i), rs_desc)
317 86526 : CALL rs_grid_set_box(grid_spme, rs=drpot(i))
318 86526 : CALL pw_multiply_with(dphi_g(i), green%p3m_charge)
319 86526 : CALL pw_transfer(dphi_g(i), rhob_r)
320 86526 : CALL pw_pool%give_back_pw(dphi_g(i))
321 115368 : CALL transfer_pw2rs(drpot(i), rhob_r)
322 : END DO
323 :
324 28842 : CALL pw_pool%give_back_pw(rhob_r)
325 28842 : DEALLOCATE (rhob_r)
326 : !----------------- FORCE CALCULATION ----------------
327 : ! initialize the forces
328 25565746 : fg_coulomb = 0.0_dp
329 : ! Particles
330 28842 : ipart = 0
331 : DO
332 3582406 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
333 3582406 : IF (p1 == 0) EXIT
334 :
335 3553564 : do_shell = (particle_set(p1)%shell_index /= 0)
336 3553564 : IF (do_shell) CYCLE
337 : ! calculate function on small boxes
338 : CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
339 3134675 : is_shell=.FALSE., unit_charge=.FALSE., charges=charges)
340 :
341 : ! add boxes to real space grid (big box)
342 3134675 : CALL dg_sum_patch_force_3d(drpot, rhos, center(:, p1), fat)
343 3134675 : fg_coulomb(1, p1) = fg_coulomb(1, p1) - fat(1)*dvols
344 3134675 : fg_coulomb(2, p1) = fg_coulomb(2, p1) - fat(2)*dvols
345 3553564 : fg_coulomb(3, p1) = fg_coulomb(3, p1) - fat(3)*dvols
346 : END DO
347 : ! Shell-Model
348 28842 : IF (PRESENT(shell_particle_set) .AND. (PRESENT(core_particle_set))) THEN
349 10396 : IF (PRESENT(fgshell_coulomb)) THEN
350 10396 : ipart = 0
351 3083348 : fgshell_coulomb = 0.0_dp
352 429285 : DO
353 : ! calculate function on small boxes
354 : CALL set_list(shell_particle_set, nshell, shell_center, p1_shell, &
355 429285 : rden, ipart)
356 429285 : IF (p1_shell == 0) EXIT
357 :
358 : CALL get_patch(shell_particle_set, shell_delta, green, &
359 418889 : p1_shell, rhos, is_core=.FALSE., is_shell=.TRUE., unit_charge=.FALSE.)
360 :
361 : ! add boxes to real space grid (big box)
362 418889 : CALL dg_sum_patch_force_3d(drpot, rhos, shell_center(:, p1_shell), fat)
363 418889 : fgshell_coulomb(1, p1_shell) = fgshell_coulomb(1, p1_shell) - fat(1)*dvols
364 418889 : fgshell_coulomb(2, p1_shell) = fgshell_coulomb(2, p1_shell) - fat(2)*dvols
365 418889 : fgshell_coulomb(3, p1_shell) = fgshell_coulomb(3, p1_shell) - fat(3)*dvols
366 :
367 : END DO
368 : END IF
369 10396 : IF (PRESENT(fgcore_coulomb)) THEN
370 10396 : ipart = 0
371 3083348 : fgcore_coulomb = 0.0_dp
372 429285 : DO
373 : ! calculate function on small boxes
374 : CALL set_list(core_particle_set, nshell, core_center, p1_shell, &
375 429285 : rden, ipart)
376 429285 : IF (p1_shell == 0) EXIT
377 :
378 : CALL get_patch(core_particle_set, core_delta, green, &
379 418889 : p1_shell, rhos, is_core=.TRUE., is_shell=.FALSE., unit_charge=.FALSE.)
380 :
381 : ! add boxes to real space grid (big box)
382 418889 : CALL dg_sum_patch_force_3d(drpot, rhos, core_center(:, p1_shell), fat)
383 418889 : fgcore_coulomb(1, p1_shell) = fgcore_coulomb(1, p1_shell) - fat(1)*dvols
384 418889 : fgcore_coulomb(2, p1_shell) = fgcore_coulomb(2, p1_shell) - fat(2)*dvols
385 418889 : 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 28842 : CALL rs_grid_release(rden)
394 28842 : DEALLOCATE (rden)
395 115368 : DO i = 1, 3
396 115368 : CALL rs_grid_release(drpot(i))
397 : END DO
398 :
399 28842 : DEALLOCATE (rhos)
400 28842 : DEALLOCATE (center, delta)
401 28842 : IF (ALLOCATED(shell_center)) THEN
402 10396 : DEALLOCATE (shell_center, shell_delta)
403 : END IF
404 28842 : IF (ALLOCATED(core_center)) THEN
405 10396 : DEALLOCATE (core_center, core_delta)
406 : END IF
407 28842 : CALL timestop(handle)
408 :
409 144210 : 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 16970 : SUBROUTINE spme_potential(ewald_env, ewald_pw, box, particle_set_a, charges_a, &
426 16970 : 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 16970 : NULLIFY (grid_spme, green, poisson_env, phi_g, rhob_g, rhob_r, pw_pool, &
456 16970 : rden, rpot)
457 16970 : CALL timeset(routineN, handle)
458 16970 : CALL cite_reference(Essmann1995)
459 :
460 : !-------------- INITIALISATION ---------------------
461 16970 : CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group)
462 16970 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, poisson_env=poisson_env)
463 16970 : CALL pw_poisson_rebuild(poisson_env)
464 16970 : green => poisson_env%green_fft
465 16970 : grid_spme => pw_pool%pw_grid
466 :
467 16970 : npart_a = SIZE(particle_set_a)
468 16970 : npart_b = SIZE(particle_set_b)
469 :
470 16970 : CALL get_pw_grid_info(grid_spme, npts=npts, dvol=dvols)
471 :
472 16970 : n = o_spline
473 84850 : ALLOCATE (rhos(n, n, n))
474 356370 : ALLOCATE (rden)
475 16970 : CALL rs_grid_create(rden, rs_desc)
476 16970 : CALL rs_grid_set_box(grid_spme, rs=rden)
477 16970 : CALL rs_grid_zero(rden)
478 :
479 84850 : ALLOCATE (center(3, npart_a), delta(3, npart_a))
480 16970 : CALL get_center(particle_set_a, box, center, delta, npts, n)
481 :
482 : !-------------- DENSITY CALCULATION ----------------
483 16970 : ipart = 0
484 : ! Particles
485 52522 : DO
486 69492 : CALL set_list(particle_set_a, npart_a, center, p1, rden, ipart)
487 69492 : IF (p1 == 0) EXIT
488 :
489 : ! calculate function on small boxes
490 52522 : CALL get_patch(particle_set_a, delta, green, p1, rhos, charges_a)
491 :
492 : ! add boxes to real space grid (big box)
493 69492 : CALL dg_sum_patch(rden, rhos, center(:, p1))
494 : END DO
495 16970 : DEALLOCATE (center, delta)
496 : !----------- END OF DENSITY CALCULATION -------------
497 :
498 16970 : NULLIFY (rhob_r)
499 16970 : ALLOCATE (rhob_r)
500 16970 : CALL pw_pool%create_pw(rhob_r)
501 16970 : CALL transfer_rs2pw(rden, rhob_r)
502 : ! transform density to G space and add charge function
503 16970 : NULLIFY (rhob_g)
504 16970 : ALLOCATE (rhob_g)
505 16970 : CALL pw_pool%create_pw(rhob_g)
506 16970 : CALL pw_transfer(rhob_r, rhob_g)
507 : ! update charge function
508 16970 : CALL pw_multiply_with(rhob_g, green%p3m_charge)
509 :
510 : !-------------- ELECTROSTATIC CALCULATION -----------
511 : ! allocate intermediate arrays
512 16970 : NULLIFY (phi_g)
513 16970 : ALLOCATE (phi_g)
514 16970 : CALL pw_pool%create_pw(phi_g)
515 16970 : CALL pw_poisson_solve(poisson_env, density=rhob_g, vhartree=phi_g)
516 : !---------- END OF ELECTROSTATIC CALCULATION --------
517 :
518 : !-------------- POTENTAIL AT ATOMIC POSITIONS -------
519 84850 : ALLOCATE (center(3, npart_b), delta(3, npart_b))
520 16970 : CALL get_center(particle_set_b, box, center, delta, npts, n)
521 16970 : rpot => rden
522 16970 : CALL pw_multiply_with(phi_g, green%p3m_charge)
523 16970 : CALL pw_transfer(phi_g, rhob_r)
524 16970 : CALL transfer_pw2rs(rpot, rhob_r)
525 16970 : ipart = 0
526 52522 : DO
527 69492 : CALL set_list(particle_set_b, npart_b, center, p1, rden, ipart)
528 69492 : IF (p1 == 0) EXIT
529 : ! calculate function on small boxes
530 : CALL get_patch(particle_set_b, delta, green, p1, rhos, &
531 52522 : is_core=.FALSE., is_shell=.FALSE., unit_charge=.TRUE.)
532 : ! integrate box and potential
533 52522 : CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fat1)
534 52522 : potential(p1) = potential(p1) + fat1*dvols
535 : END DO
536 :
537 : !------------------CLEANING UP ----------------------
538 16970 : CALL pw_pool%give_back_pw(phi_g)
539 16970 : CALL pw_pool%give_back_pw(rhob_g)
540 16970 : CALL pw_pool%give_back_pw(rhob_r)
541 16970 : DEALLOCATE (phi_g, rhob_g, rhob_r)
542 16970 : CALL rs_grid_release(rden)
543 16970 : DEALLOCATE (rden)
544 :
545 16970 : DEALLOCATE (rhos)
546 16970 : DEALLOCATE (center, delta)
547 16970 : CALL timestop(handle)
548 :
549 50910 : 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 192 : SUBROUTINE spme_forces(ewald_env, ewald_pw, box, particle_set_a, charges_a, &
567 192 : 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 192 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: center
583 : INTEGER, DIMENSION(3) :: npts
584 : REAL(KIND=dp) :: alpha, dvols
585 192 : 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 768 : 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 6912 : TYPE(realspace_grid_type), DIMENSION(3) :: drpot
598 : TYPE(realspace_grid_type), POINTER :: rden
599 :
600 192 : NULLIFY (grid_spme, green, poisson_env, phi_g, rhob_g, rhob_r, &
601 192 : pw_pool, rden)
602 192 : CALL timeset(routineN, handle)
603 192 : CALL cite_reference(Essmann1995)
604 :
605 : !-------------- INITIALISATION ---------------------
606 192 : CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group)
607 192 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, poisson_env=poisson_env)
608 192 : CALL pw_poisson_rebuild(poisson_env)
609 192 : green => poisson_env%green_fft
610 192 : grid_spme => pw_pool%pw_grid
611 :
612 192 : npart_a = SIZE(particle_set_a)
613 192 : npart_b = SIZE(particle_set_b)
614 :
615 192 : CALL get_pw_grid_info(grid_spme, npts=npts, dvol=dvols)
616 :
617 192 : n = o_spline
618 960 : ALLOCATE (rhos(n, n, n))
619 4032 : ALLOCATE (rden)
620 192 : CALL rs_grid_create(rden, rs_desc)
621 192 : CALL rs_grid_set_box(grid_spme, rs=rden)
622 192 : CALL rs_grid_zero(rden)
623 :
624 960 : ALLOCATE (center(3, npart_a), delta(3, npart_a))
625 192 : CALL get_center(particle_set_a, box, center, delta, npts, n)
626 :
627 : !-------------- DENSITY CALCULATION ----------------
628 192 : ipart = 0
629 : ! Particles
630 452 : DO
631 644 : CALL set_list(particle_set_a, npart_a, center, p1, rden, ipart)
632 644 : IF (p1 == 0) EXIT
633 :
634 : ! calculate function on small boxes
635 452 : CALL get_patch(particle_set_a, delta, green, p1, rhos, charges_a)
636 :
637 : ! add boxes to real space grid (big box)
638 644 : CALL dg_sum_patch(rden, rhos, center(:, p1))
639 : END DO
640 192 : DEALLOCATE (center, delta)
641 : !----------- END OF DENSITY CALCULATION -------------
642 :
643 192 : NULLIFY (rhob_r)
644 192 : ALLOCATE (rhob_r)
645 192 : CALL pw_pool%create_pw(rhob_r)
646 192 : CALL transfer_rs2pw(rden, rhob_r)
647 : ! transform density to G space and add charge function
648 192 : NULLIFY (rhob_g)
649 192 : ALLOCATE (rhob_g)
650 192 : CALL pw_pool%create_pw(rhob_g)
651 192 : CALL pw_transfer(rhob_r, rhob_g)
652 : ! update charge function
653 192 : CALL pw_multiply_with(rhob_g, green%p3m_charge)
654 :
655 : !-------------- ELECTROSTATIC CALCULATION -----------
656 : ! allocate intermediate arrays
657 768 : DO i = 1, 3
658 768 : CALL pw_pool%create_pw(dphi_g(i))
659 : END DO
660 192 : NULLIFY (phi_g)
661 192 : ALLOCATE (phi_g)
662 192 : CALL pw_pool%create_pw(phi_g)
663 : CALL pw_poisson_solve(poisson_env, density=rhob_g, vhartree=phi_g, &
664 192 : 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 768 : DO i = 1, 3
669 576 : CALL rs_grid_create(drpot(i), rs_desc)
670 576 : CALL rs_grid_set_box(grid_spme, rs=drpot(i))
671 576 : CALL pw_multiply_with(dphi_g(i), green%p3m_charge)
672 576 : CALL pw_transfer(dphi_g(i), rhob_r)
673 576 : CALL pw_pool%give_back_pw(dphi_g(i))
674 768 : CALL transfer_pw2rs(drpot(i), rhob_r)
675 : END DO
676 : !----------------- FORCE CALCULATION ----------------
677 960 : ALLOCATE (center(3, npart_b), delta(3, npart_b))
678 192 : CALL get_center(particle_set_b, box, center, delta, npts, n)
679 192 : ipart = 0
680 452 : DO
681 644 : CALL set_list(particle_set_b, npart_b, center, p1, rden, ipart)
682 644 : IF (p1 == 0) EXIT
683 : ! calculate function on small boxes
684 : CALL get_patch(particle_set_b, delta, green, p1, rhos, &
685 452 : is_core=.FALSE., is_shell=.FALSE., unit_charge=.FALSE., charges=charges_b)
686 : ! add boxes to real space grid (big box)
687 452 : CALL dg_sum_patch_force_3d(drpot, rhos, center(:, p1), fat)
688 452 : forces_b(1, p1) = forces_b(1, p1) - fat(1)*dvols
689 452 : forces_b(2, p1) = forces_b(2, p1) - fat(2)*dvols
690 644 : forces_b(3, p1) = forces_b(3, p1) - fat(3)*dvols
691 : END DO
692 : !------------------CLEANING UP ----------------------
693 768 : DO i = 1, 3
694 768 : CALL rs_grid_release(drpot(i))
695 : END DO
696 192 : CALL pw_pool%give_back_pw(phi_g)
697 192 : CALL pw_pool%give_back_pw(rhob_g)
698 192 : CALL pw_pool%give_back_pw(rhob_r)
699 192 : DEALLOCATE (phi_g, rhob_g, rhob_r)
700 192 : CALL rs_grid_release(rden)
701 192 : DEALLOCATE (rden)
702 :
703 192 : DEALLOCATE (rhos)
704 192 : DEALLOCATE (center, delta)
705 192 : CALL timestop(handle)
706 :
707 2880 : 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 72 : 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 72 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: center
732 : INTEGER, DIMENSION(3) :: npts
733 : REAL(KIND=dp) :: alpha, dvols, ffa, vgc
734 72 : 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 288 : 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 2808 : TYPE(realspace_grid_type) :: rden, rpot
748 :
749 72 : CALL timeset(routineN, handle)
750 : !-------------- INITIALISATION ---------------------
751 72 : virial = 0.0_dp
752 : CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group, &
753 72 : para_env=para_env)
754 72 : NULLIFY (green, poisson_env, pw_pool)
755 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, &
756 72 : poisson_env=poisson_env)
757 72 : CALL pw_poisson_rebuild(poisson_env)
758 72 : green => poisson_env%green_fft
759 72 : grid_spme => pw_pool%pw_grid
760 :
761 72 : CALL get_pw_grid_info(grid_spme, dvol=dvols, npts=npts)
762 :
763 72 : npart = SIZE(particle_set)
764 :
765 72 : n = o_spline
766 360 : ALLOCATE (rhos(n, n, n))
767 :
768 72 : CALL rs_grid_create(rden, rs_desc)
769 72 : CALL rs_grid_set_box(grid_spme, rs=rden)
770 72 : CALL rs_grid_zero(rden)
771 :
772 360 : ALLOCATE (center(3, npart), delta(3, npart))
773 72 : CALL get_center(particle_set, box, center, delta, npts, n)
774 :
775 : !-------------- DENSITY CALCULATION ----------------
776 72 : ipart = 0
777 240 : DO
778 312 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
779 312 : 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 240 : is_shell=.FALSE., unit_charge=.TRUE.)
784 20400 : rhos(:, :, :) = rhos(:, :, :)*mcharge(p1)
785 :
786 : ! add boxes to real space grid (big box)
787 312 : CALL dg_sum_patch(rden, rhos, center(:, p1))
788 : END DO
789 :
790 72 : NULLIFY (rhob_r)
791 72 : ALLOCATE (rhob_r)
792 72 : CALL pw_pool%create_pw(rhob_r)
793 :
794 72 : CALL transfer_rs2pw(rden, rhob_r)
795 :
796 : ! transform density to G space and add charge function
797 72 : NULLIFY (rhob_g)
798 72 : ALLOCATE (rhob_g)
799 72 : CALL pw_pool%create_pw(rhob_g)
800 72 : CALL pw_transfer(rhob_r, rhob_g)
801 : ! update charge function
802 72 : CALL pw_multiply_with(rhob_g, green%p3m_charge)
803 :
804 : !-------------- ELECTROSTATIC CALCULATION -----------
805 :
806 : ! allocate intermediate arrays
807 288 : DO i = 1, 3
808 288 : CALL pw_pool%create_pw(dphi_g(i))
809 : END DO
810 72 : NULLIFY (phi_g)
811 72 : ALLOCATE (phi_g)
812 72 : CALL pw_pool%create_pw(phi_g)
813 72 : CALL pw_poisson_solve(poisson_env, rhob_g, vgc, phi_g, dphi_g, h_stress=h_stress)
814 :
815 72 : CALL rs_grid_create(rpot, rs_desc)
816 72 : CALL rs_grid_set_box(grid_spme, rs=rpot)
817 :
818 72 : CALL pw_pool%give_back_pw(rhob_g)
819 72 : DEALLOCATE (rhob_g)
820 :
821 72 : CALL rs_grid_zero(rpot)
822 72 : CALL pw_multiply_with(phi_g, green%p3m_charge)
823 72 : CALL pw_transfer(phi_g, rhob_r)
824 72 : CALL pw_pool%give_back_pw(phi_g)
825 72 : DEALLOCATE (phi_g)
826 72 : CALL transfer_pw2rs(rpot, rhob_r)
827 :
828 : !---------- END OF ELECTROSTATIC CALCULATION --------
829 :
830 : !------------- STRESS TENSOR CALCULATION ------------
831 :
832 288 : DO i = 1, 3
833 720 : DO j = i, 3
834 432 : f_stress(i, j) = pw_integral_a2b(dphi_g(i), dphi_g(j))
835 648 : f_stress(j, i) = f_stress(i, j)
836 : END DO
837 : END DO
838 72 : ffa = (1.0_dp/fourpi)*(0.5_dp/alpha)**2
839 936 : virial = virial - (ffa*f_stress - h_stress)/REAL(para_env%num_pe, dp)
840 :
841 : !--------END OF STRESS TENSOR CALCULATION -----------
842 :
843 288 : DO i = 1, 3
844 288 : CALL pw_pool%give_back_pw(dphi_g(i))
845 : END DO
846 72 : CALL pw_pool%give_back_pw(rhob_r)
847 72 : DEALLOCATE (rhob_r)
848 :
849 72 : CALL rs_grid_release(rden)
850 72 : CALL rs_grid_release(rpot)
851 72 : DEALLOCATE (rhos)
852 72 : DEALLOCATE (center, delta)
853 :
854 72 : CALL timestop(handle)
855 :
856 216 : 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 8322213 : 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 8322213 : NULLIFY (shell)
891 8322213 : use_charge_array = .FALSE.
892 8322213 : IF (PRESENT(charges)) use_charge_array = ASSOCIATED(charges)
893 8322213 : IF (is_core .AND. is_shell) THEN
894 0 : CPABORT("Shell-model: cannot be core and shell simultaneously")
895 : END IF
896 :
897 8322213 : nbox = SIZE(rhos, 1)
898 33288852 : r = part(p)%r
899 8322213 : q = 1.0_dp
900 8322213 : IF (.NOT. unit_charge) THEN
901 7948541 : IF (is_core) THEN
902 838858 : CALL get_atomic_kind(atomic_kind=part(p)%atomic_kind, shell=shell)
903 838858 : q = shell%charge_core
904 7109683 : ELSE IF (is_shell) THEN
905 838858 : CALL get_atomic_kind(atomic_kind=part(p)%atomic_kind, shell=shell)
906 838858 : q = shell%charge_shell
907 : ELSE
908 6270825 : CALL get_atomic_kind(atomic_kind=part(p)%atomic_kind, qeff=q)
909 : END IF
910 7948541 : IF (use_charge_array) q = charges(p)
911 : END IF
912 8322213 : CALL spme_get_patch(rhos, nbox, delta(:, p), q, green%p3m_coeff)
913 :
914 8322213 : 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 52974 : 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 52974 : CPASSERT(ASSOCIATED(charges))
942 52974 : nbox = SIZE(rhos, 1)
943 211896 : r = part(p)%r
944 52974 : q = charges(p)
945 52974 : CALL spme_get_patch(rhos, nbox, delta(:, p), q, green%p3m_coeff)
946 :
947 52974 : 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 8375187 : 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 8375187 : 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 8375187 : deltal(1, 0) = 1.0_dp
984 8375187 : deltal(2, 0) = 1.0_dp
985 8375187 : deltal(3, 0) = 1.0_dp
986 44089014 : DO l = 1, n - 1
987 35713827 : deltal(1, l) = deltal(1, l - 1)*delta(1)
988 35713827 : deltal(2, l) = deltal(2, l - 1)*delta(2)
989 44089014 : deltal(3, l) = deltal(3, l - 1)*delta(3)
990 : END DO
991 :
992 8375187 : w_assign = 0.0_dp
993 52464201 : DO j = -(n - 1), n - 1, 2
994 292441277 : DO l = 0, n - 1
995 239977076 : w_assign(1, j) = w_assign(1, j) + coeff(j, l)*deltal(1, l)
996 239977076 : w_assign(2, j) = w_assign(2, j) + coeff(j, l)*deltal(2, l)
997 284066090 : w_assign(3, j) = w_assign(3, j) + coeff(j, l)*deltal(3, l)
998 : END DO
999 : END DO
1000 52464201 : DO i = 1, n
1001 44089014 : j = n + 1 - 2*i
1002 44089014 : f_assign(1, i) = w_assign(1, j)
1003 44089014 : f_assign(2, i) = w_assign(2, j)
1004 52464201 : f_assign(3, i) = w_assign(3, j)
1005 : END DO
1006 :
1007 52464201 : DO i3 = 1, n
1008 44089014 : r3 = q*f_assign(3, i3)
1009 292441277 : DO i2 = 1, n
1010 239977076 : r2 = r3*f_assign(2, i2)
1011 1626641474 : DO i1 = 1, n
1012 1582552460 : rhos(i1, i2, i3) = r2*f_assign(1, i1)
1013 : END DO
1014 : END DO
1015 : END DO
1016 :
1017 8375187 : END SUBROUTINE spme_get_patch
1018 :
1019 : END MODULE spme
1020 :
|
