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