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