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