(git:374b731)
Loading...
Searching...
No Matches
spme.F
Go to the documentation of this file.
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! **************************************************************************************************
14MODULE spme
15
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,&
26 USE ewald_pw_types, ONLY: ewald_pw_get,&
28 USE kinds, ONLY: dp
29 USE mathconstants, ONLY: fourpi
32 USE pme_tools, ONLY: get_center,&
35 USE pw_grids, ONLY: get_pw_grid_info
36 USE pw_methods, ONLY: pw_copy,&
37 pw_derive,&
46 USE pw_types, ONLY: pw_c1d_gs_type,&
57#include "./base/base_uses.f90"
58
59 IMPLICIT NONE
60
61 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'spme'
62
63 PRIVATE
65
66 INTERFACE get_patch
67 MODULE PROCEDURE get_patch_a, get_patch_b
68 END INTERFACE
69
70CONTAINS
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 SUBROUTINE spme_evaluate(ewald_env, ewald_pw, box, particle_set, &
93 fg_coulomb, vg_coulomb, pv_g, shell_particle_set, core_particle_set, &
94 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 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 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 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 TYPE(realspace_grid_type), DIMENSION(3) :: drpot
133 TYPE(realspace_grid_type), POINTER :: rden, rpot
134
135 NULLIFY (grid_spme, green, poisson_env, phi_g, rhob_g, rhob_r, pw_pool, rden, rpot)
136 CALL timeset(routinen, handle)
137 CALL cite_reference(essmann1995)
138
139 !-------------- INITIALISATION ---------------------
140 CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, &
141 group=group)
142 CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, &
143 poisson_env=poisson_env)
144 CALL pw_poisson_rebuild(poisson_env)
145 green => poisson_env%green_fft
146 grid_spme => pw_pool%pw_grid
147
148 npart = SIZE(particle_set)
149
150 CALL get_pw_grid_info(grid_spme, npts=npts, dvol=dvols)
151
152 n = o_spline
153 ALLOCATE (rhos(n, n, n))
154 ALLOCATE (rden)
155 CALL rs_grid_create(rden, rs_desc)
156 CALL rs_grid_set_box(grid_spme, rs=rden)
157 CALL rs_grid_zero(rden)
158
159 ALLOCATE (center(3, npart), delta(3, npart))
160 CALL get_center(particle_set, box, center, delta, npts, n)
161
162 IF (PRESENT(shell_particle_set) .AND. (PRESENT(core_particle_set))) THEN
163 cpassert(ASSOCIATED(shell_particle_set))
164 cpassert(ASSOCIATED(core_particle_set))
165 nshell = SIZE(shell_particle_set)
166 ncore = SIZE(core_particle_set)
167 cpassert(nshell == ncore)
168 ALLOCATE (shell_center(3, nshell), shell_delta(3, nshell))
169 CALL get_center(shell_particle_set, box, shell_center, shell_delta, npts, n)
170 ALLOCATE (core_center(3, nshell), core_delta(3, nshell))
171 CALL get_center(core_particle_set, box, core_center, core_delta, npts, n)
172 END IF
173
174 !-------------- DENSITY CALCULATION ----------------
175 ipart = 0
176 ! Particles
177 DO
178 CALL set_list(particle_set, npart, center, p1, rden, ipart)
179 IF (p1 == 0) EXIT
180
181 do_shell = (particle_set(p1)%shell_index /= 0)
182 IF (do_shell) cycle
183 ! calculate function on small boxes
184 CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.false., &
185 is_shell=.false., unit_charge=.false., charges=charges)
186
187 ! add boxes to real space grid (big box)
188 CALL dg_sum_patch(rden, rhos, center(:, p1))
189 END DO
190 ! Shell-Model
191 IF (PRESENT(shell_particle_set) .AND. PRESENT(core_particle_set)) THEN
192 ipart = 0
193 DO
194 ! calculate function on small boxes
195 CALL set_list(shell_particle_set, nshell, shell_center, p1_shell, &
196 rden, ipart)
197 IF (p1_shell == 0) EXIT
198 CALL get_patch(shell_particle_set, shell_delta, green, p1_shell, rhos, &
199 is_core=.false., is_shell=.true., unit_charge=.false.)
200
201 ! add boxes to real space grid (big box)
202 CALL dg_sum_patch(rden, rhos, shell_center(:, p1_shell))
203 END DO
204 ipart = 0
205 DO
206 ! calculate function on small boxes
207 CALL set_list(core_particle_set, nshell, core_center, p1_shell, &
208 rden, ipart)
209 IF (p1_shell == 0) EXIT
210 CALL get_patch(core_particle_set, core_delta, green, p1_shell, rhos, &
211 is_core=.true., is_shell=.false., unit_charge=.false.)
212
213 ! add boxes to real space grid (big box)
214 CALL dg_sum_patch(rden, rhos, core_center(:, p1_shell))
215 END DO
216 END IF
217 !----------- END OF DENSITY CALCULATION -------------
218
219 NULLIFY (rhob_r)
220 ALLOCATE (rhob_r)
221 CALL pw_pool%create_pw(rhob_r)
222 CALL transfer_rs2pw(rden, rhob_r)
223 ! transform density to G space and add charge function
224 NULLIFY (rhob_g)
225 ALLOCATE (rhob_g)
226 CALL pw_pool%create_pw(rhob_g)
227 CALL pw_transfer(rhob_r, rhob_g)
228 ! update charge function
229 CALL pw_multiply_with(rhob_g, green%p3m_charge)
230
231 !-------------- ELECTROSTATIC CALCULATION -----------
232 ! allocate intermediate arrays
233 DO i = 1, 3
234 CALL pw_pool%create_pw(dphi_g(i))
235 END DO
236 NULLIFY (phi_g)
237 ALLOCATE (phi_g)
238 CALL pw_pool%create_pw(phi_g)
239 CALL pw_poisson_solve(poisson_env, rhob_g, vg_coulomb, phi_g, dphi_g, &
240 h_stress)
241 !---------- END OF ELECTROSTATIC CALCULATION --------
242
243 ! Atomic Energy and Stress
244 IF (atprop%energy .OR. atprop%stress) THEN
245 ALLOCATE (rpot)
246 CALL rs_grid_create(rpot, rs_desc)
247 CALL rs_grid_set_box(grid_spme, rs=rpot)
248 CALL pw_multiply_with(phi_g, green%p3m_charge)
249 CALL pw_transfer(phi_g, rhob_r)
250 CALL transfer_pw2rs(rpot, rhob_r)
251 ipart = 0
252 DO
253 CALL set_list(particle_set, npart, center, p1, rden, ipart)
254 IF (p1 == 0) EXIT
255 do_shell = (particle_set(p1)%shell_index /= 0)
256 IF (do_shell) cycle
257 ! calculate function on small boxes
258 CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.false., &
259 is_shell=.false., unit_charge=.false., charges=charges)
260 ! integrate box and potential
261 CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fat1)
262 IF (atprop%energy) THEN
263 atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols
264 END IF
265 IF (atprop%stress) THEN
266 atprop%atstress(1, 1, p1) = atprop%atstress(1, 1, p1) + 0.5_dp*fat1*dvols
267 atprop%atstress(2, 2, p1) = atprop%atstress(2, 2, p1) + 0.5_dp*fat1*dvols
268 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 IF (PRESENT(shell_particle_set) .AND. PRESENT(core_particle_set)) THEN
273 ipart = 0
274 DO
275 CALL set_list(shell_particle_set, nshell, shell_center, p1_shell, rden, ipart)
276 IF (p1_shell == 0) EXIT
277 CALL get_patch(shell_particle_set, shell_delta, green, p1_shell, rhos, &
278 is_core=.false., is_shell=.true., unit_charge=.false.)
279 CALL dg_sum_patch_force_1d(rpot, rhos, shell_center(:, p1_shell), fat1)
280 p1 = shell_particle_set(p1_shell)%atom_index
281 IF (atprop%energy) THEN
282 atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols
283 END IF
284 END DO
285 ipart = 0
286 DO
287 CALL set_list(core_particle_set, nshell, core_center, p1_shell, rden, ipart)
288 IF (p1_shell == 0) EXIT
289 CALL get_patch(core_particle_set, core_delta, green, p1_shell, rhos, &
290 is_core=.true., is_shell=.false., unit_charge=.false.)
291 CALL dg_sum_patch_force_1d(rpot, rhos, core_center(:, p1_shell), fat1)
292 p1 = core_particle_set(p1_shell)%atom_index
293 IF (atprop%energy) THEN
294 atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols
295 END IF
296 END DO
297 END IF
298 IF (atprop%stress) THEN
299 ffa = (0.5_dp/alpha)**2
300 ffb = 1.0_dp/fourpi
301 DO i = 1, 3
302 DO ig = grid_spme%first_gne0, grid_spme%ngpts_cut_local
303 phi_g%array(ig) = ffb*dphi_g(i)%array(ig)*(ffa*grid_spme%gsq(ig) + 1.0_dp)
304 phi_g%array(ig) = phi_g%array(ig)*poisson_env%green_fft%influence_fn%array(ig)
305 END DO
306 IF (grid_spme%have_g0) phi_g%array(1) = 0.0_dp
307 DO j = 1, i
308 nd = 0
309 nd(j) = 1
310 CALL pw_copy(phi_g, rhob_g)
311 CALL pw_derive(rhob_g, nd)
312 CALL pw_multiply_with(rhob_g, green%p3m_charge)
313 CALL pw_transfer(rhob_g, rhob_r)
314 CALL transfer_pw2rs(rpot, rhob_r)
315
316 ipart = 0
317 DO
318 CALL set_list(particle_set, npart, center, p1, rden, ipart)
319 IF (p1 == 0) EXIT
320 ! calculate function on small boxes
321 CALL get_patch(particle_set, delta, green, p1, rhos, &
322 is_core=.false., is_shell=.false., unit_charge=.false., charges=charges)
323 ! integrate box and potential
324 CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fat1)
325 atprop%atstress(i, j, p1) = atprop%atstress(i, j, p1) + fat1*dvols
326 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 CALL rs_grid_release(rpot)
333 DEALLOCATE (rpot)
334 END IF
335
336 CALL pw_pool%give_back_pw(phi_g)
337 CALL pw_pool%give_back_pw(rhob_g)
338 DEALLOCATE (phi_g, rhob_g)
339
340 !------------- STRESS TENSOR CALCULATION ------------
341 IF (use_virial) THEN
342 DO i = 1, 3
343 DO j = i, 3
344 f_stress(i, j) = pw_integral_a2b(dphi_g(i), dphi_g(j))
345 f_stress(j, i) = f_stress(i, j)
346 END DO
347 END DO
348 ffa = (1.0_dp/fourpi)*(0.5_dp/alpha)**2
349 f_stress = -ffa*f_stress
350 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 DO i = 1, 3
356 CALL rs_grid_create(drpot(i), rs_desc)
357 CALL rs_grid_set_box(grid_spme, rs=drpot(i))
358 CALL pw_multiply_with(dphi_g(i), green%p3m_charge)
359 CALL pw_transfer(dphi_g(i), rhob_r)
360 CALL pw_pool%give_back_pw(dphi_g(i))
361 CALL transfer_pw2rs(drpot(i), rhob_r)
362 END DO
363
364 CALL pw_pool%give_back_pw(rhob_r)
365 DEALLOCATE (rhob_r)
366 !----------------- FORCE CALCULATION ----------------
367 ! initialize the forces
368 fg_coulomb = 0.0_dp
369 ! Particles
370 ipart = 0
371 DO
372 CALL set_list(particle_set, npart, center, p1, rden, ipart)
373 IF (p1 == 0) EXIT
374
375 do_shell = (particle_set(p1)%shell_index /= 0)
376 IF (do_shell) cycle
377 ! calculate function on small boxes
378 CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.false., &
379 is_shell=.false., unit_charge=.false., charges=charges)
380
381 ! add boxes to real space grid (big box)
382 CALL dg_sum_patch_force_3d(drpot, rhos, center(:, p1), fat)
383 fg_coulomb(1, p1) = fg_coulomb(1, p1) - fat(1)*dvols
384 fg_coulomb(2, p1) = fg_coulomb(2, p1) - fat(2)*dvols
385 fg_coulomb(3, p1) = fg_coulomb(3, p1) - fat(3)*dvols
386 END DO
387 ! Shell-Model
388 IF (PRESENT(shell_particle_set) .AND. (PRESENT(core_particle_set))) THEN
389 IF (PRESENT(fgshell_coulomb)) THEN
390 ipart = 0
391 fgshell_coulomb = 0.0_dp
392 DO
393 ! calculate function on small boxes
394 CALL set_list(shell_particle_set, nshell, shell_center, p1_shell, &
395 rden, ipart)
396 IF (p1_shell == 0) EXIT
397
398 CALL get_patch(shell_particle_set, shell_delta, green, &
399 p1_shell, rhos, is_core=.false., is_shell=.true., unit_charge=.false.)
400
401 ! add boxes to real space grid (big box)
402 CALL dg_sum_patch_force_3d(drpot, rhos, shell_center(:, p1_shell), fat)
403 fgshell_coulomb(1, p1_shell) = fgshell_coulomb(1, p1_shell) - fat(1)*dvols
404 fgshell_coulomb(2, p1_shell) = fgshell_coulomb(2, p1_shell) - fat(2)*dvols
405 fgshell_coulomb(3, p1_shell) = fgshell_coulomb(3, p1_shell) - fat(3)*dvols
406
407 END DO
408 END IF
409 IF (PRESENT(fgcore_coulomb)) THEN
410 ipart = 0
411 fgcore_coulomb = 0.0_dp
412 DO
413 ! calculate function on small boxes
414 CALL set_list(core_particle_set, nshell, core_center, p1_shell, &
415 rden, ipart)
416 IF (p1_shell == 0) EXIT
417
418 CALL get_patch(core_particle_set, core_delta, green, &
419 p1_shell, rhos, is_core=.true., is_shell=.false., unit_charge=.false.)
420
421 ! add boxes to real space grid (big box)
422 CALL dg_sum_patch_force_3d(drpot, rhos, core_center(:, p1_shell), fat)
423 fgcore_coulomb(1, p1_shell) = fgcore_coulomb(1, p1_shell) - fat(1)*dvols
424 fgcore_coulomb(2, p1_shell) = fgcore_coulomb(2, p1_shell) - fat(2)*dvols
425 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 CALL rs_grid_release(rden)
434 DEALLOCATE (rden)
435 DO i = 1, 3
436 CALL rs_grid_release(drpot(i))
437 END DO
438
439 DEALLOCATE (rhos)
440 DEALLOCATE (center, delta)
441 IF (ALLOCATED(shell_center)) THEN
442 DEALLOCATE (shell_center, shell_delta)
443 END IF
444 IF (ALLOCATED(core_center)) THEN
445 DEALLOCATE (core_center, core_delta)
446 END IF
447 CALL timestop(handle)
448
449 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 SUBROUTINE spme_potential(ewald_env, ewald_pw, box, particle_set_a, charges_a, &
466 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 NULLIFY (grid_spme, green, poisson_env, phi_g, rhob_g, rhob_r, pw_pool, &
496 rden, rpot)
497 CALL timeset(routinen, handle)
498 CALL cite_reference(essmann1995)
499
500 !-------------- INITIALISATION ---------------------
501 CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group)
502 CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, poisson_env=poisson_env)
503 CALL pw_poisson_rebuild(poisson_env)
504 green => poisson_env%green_fft
505 grid_spme => pw_pool%pw_grid
506
507 npart_a = SIZE(particle_set_a)
508 npart_b = SIZE(particle_set_b)
509
510 CALL get_pw_grid_info(grid_spme, npts=npts, dvol=dvols)
511
512 n = o_spline
513 ALLOCATE (rhos(n, n, n))
514 ALLOCATE (rden)
515 CALL rs_grid_create(rden, rs_desc)
516 CALL rs_grid_set_box(grid_spme, rs=rden)
517 CALL rs_grid_zero(rden)
518
519 ALLOCATE (center(3, npart_a), delta(3, npart_a))
520 CALL get_center(particle_set_a, box, center, delta, npts, n)
521
522 !-------------- DENSITY CALCULATION ----------------
523 ipart = 0
524 ! Particles
525 DO
526 CALL set_list(particle_set_a, npart_a, center, p1, rden, ipart)
527 IF (p1 == 0) EXIT
528
529 ! calculate function on small boxes
530 CALL get_patch(particle_set_a, delta, green, p1, rhos, charges_a)
531
532 ! add boxes to real space grid (big box)
533 CALL dg_sum_patch(rden, rhos, center(:, p1))
534 END DO
535 DEALLOCATE (center, delta)
536 !----------- END OF DENSITY CALCULATION -------------
537
538 NULLIFY (rhob_r)
539 ALLOCATE (rhob_r)
540 CALL pw_pool%create_pw(rhob_r)
541 CALL transfer_rs2pw(rden, rhob_r)
542 ! transform density to G space and add charge function
543 NULLIFY (rhob_g)
544 ALLOCATE (rhob_g)
545 CALL pw_pool%create_pw(rhob_g)
546 CALL pw_transfer(rhob_r, rhob_g)
547 ! update charge function
548 CALL pw_multiply_with(rhob_g, green%p3m_charge)
549
550 !-------------- ELECTROSTATIC CALCULATION -----------
551 ! allocate intermediate arrays
552 NULLIFY (phi_g)
553 ALLOCATE (phi_g)
554 CALL pw_pool%create_pw(phi_g)
555 CALL pw_poisson_solve(poisson_env, density=rhob_g, vhartree=phi_g)
556 !---------- END OF ELECTROSTATIC CALCULATION --------
557
558 !-------------- POTENTAIL AT ATOMIC POSITIONS -------
559 ALLOCATE (center(3, npart_b), delta(3, npart_b))
560 CALL get_center(particle_set_b, box, center, delta, npts, n)
561 rpot => rden
562 CALL pw_multiply_with(phi_g, green%p3m_charge)
563 CALL pw_transfer(phi_g, rhob_r)
564 CALL transfer_pw2rs(rpot, rhob_r)
565 ipart = 0
566 DO
567 CALL set_list(particle_set_b, npart_b, center, p1, rden, ipart)
568 IF (p1 == 0) EXIT
569 ! calculate function on small boxes
570 CALL get_patch(particle_set_b, delta, green, p1, rhos, &
571 is_core=.false., is_shell=.false., unit_charge=.true.)
572 ! integrate box and potential
573 CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fat1)
574 potential(p1) = potential(p1) + fat1*dvols
575 END DO
576
577 !------------------CLEANING UP ----------------------
578 CALL pw_pool%give_back_pw(phi_g)
579 CALL pw_pool%give_back_pw(rhob_g)
580 CALL pw_pool%give_back_pw(rhob_r)
581 DEALLOCATE (phi_g, rhob_g, rhob_r)
582 CALL rs_grid_release(rden)
583 DEALLOCATE (rden)
584
585 DEALLOCATE (rhos)
586 DEALLOCATE (center, delta)
587 CALL timestop(handle)
588
589 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 SUBROUTINE spme_forces(ewald_env, ewald_pw, box, particle_set_a, charges_a, &
607 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 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: center
623 INTEGER, DIMENSION(3) :: npts
624 REAL(kind=dp) :: alpha, dvols
625 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 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 TYPE(realspace_grid_type), DIMENSION(3) :: drpot
638 TYPE(realspace_grid_type), POINTER :: rden
639
640 NULLIFY (grid_spme, green, poisson_env, phi_g, rhob_g, rhob_r, &
641 pw_pool, rden)
642 CALL timeset(routinen, handle)
643 CALL cite_reference(essmann1995)
644
645 !-------------- INITIALISATION ---------------------
646 CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group)
647 CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, poisson_env=poisson_env)
648 CALL pw_poisson_rebuild(poisson_env)
649 green => poisson_env%green_fft
650 grid_spme => pw_pool%pw_grid
651
652 npart_a = SIZE(particle_set_a)
653 npart_b = SIZE(particle_set_b)
654
655 CALL get_pw_grid_info(grid_spme, npts=npts, dvol=dvols)
656
657 n = o_spline
658 ALLOCATE (rhos(n, n, n))
659 ALLOCATE (rden)
660 CALL rs_grid_create(rden, rs_desc)
661 CALL rs_grid_set_box(grid_spme, rs=rden)
662 CALL rs_grid_zero(rden)
663
664 ALLOCATE (center(3, npart_a), delta(3, npart_a))
665 CALL get_center(particle_set_a, box, center, delta, npts, n)
666
667 !-------------- DENSITY CALCULATION ----------------
668 ipart = 0
669 ! Particles
670 DO
671 CALL set_list(particle_set_a, npart_a, center, p1, rden, ipart)
672 IF (p1 == 0) EXIT
673
674 ! calculate function on small boxes
675 CALL get_patch(particle_set_a, delta, green, p1, rhos, charges_a)
676
677 ! add boxes to real space grid (big box)
678 CALL dg_sum_patch(rden, rhos, center(:, p1))
679 END DO
680 DEALLOCATE (center, delta)
681 !----------- END OF DENSITY CALCULATION -------------
682
683 NULLIFY (rhob_r)
684 ALLOCATE (rhob_r)
685 CALL pw_pool%create_pw(rhob_r)
686 CALL transfer_rs2pw(rden, rhob_r)
687 ! transform density to G space and add charge function
688 NULLIFY (rhob_g)
689 ALLOCATE (rhob_g)
690 CALL pw_pool%create_pw(rhob_g)
691 CALL pw_transfer(rhob_r, rhob_g)
692 ! update charge function
693 CALL pw_multiply_with(rhob_g, green%p3m_charge)
694
695 !-------------- ELECTROSTATIC CALCULATION -----------
696 ! allocate intermediate arrays
697 DO i = 1, 3
698 CALL pw_pool%create_pw(dphi_g(i))
699 END DO
700 NULLIFY (phi_g)
701 ALLOCATE (phi_g)
702 CALL pw_pool%create_pw(phi_g)
703 CALL pw_poisson_solve(poisson_env, density=rhob_g, vhartree=phi_g, &
704 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 DO i = 1, 3
709 CALL rs_grid_create(drpot(i), rs_desc)
710 CALL rs_grid_set_box(grid_spme, rs=drpot(i))
711 CALL pw_multiply_with(dphi_g(i), green%p3m_charge)
712 CALL pw_transfer(dphi_g(i), rhob_r)
713 CALL pw_pool%give_back_pw(dphi_g(i))
714 CALL transfer_pw2rs(drpot(i), rhob_r)
715 END DO
716 !----------------- FORCE CALCULATION ----------------
717 ALLOCATE (center(3, npart_b), delta(3, npart_b))
718 CALL get_center(particle_set_b, box, center, delta, npts, n)
719 ipart = 0
720 DO
721 CALL set_list(particle_set_b, npart_b, center, p1, rden, ipart)
722 IF (p1 == 0) EXIT
723 ! calculate function on small boxes
724 CALL get_patch(particle_set_b, delta, green, p1, rhos, &
725 is_core=.false., is_shell=.false., unit_charge=.false., charges=charges_b)
726 ! add boxes to real space grid (big box)
727 CALL dg_sum_patch_force_3d(drpot, rhos, center(:, p1), fat)
728 forces_b(1, p1) = forces_b(1, p1) - fat(1)*dvols
729 forces_b(2, p1) = forces_b(2, p1) - fat(2)*dvols
730 forces_b(3, p1) = forces_b(3, p1) - fat(3)*dvols
731 END DO
732 !------------------CLEANING UP ----------------------
733 DO i = 1, 3
734 CALL rs_grid_release(drpot(i))
735 END DO
736 CALL pw_pool%give_back_pw(phi_g)
737 CALL pw_pool%give_back_pw(rhob_g)
738 CALL pw_pool%give_back_pw(rhob_r)
739 DEALLOCATE (phi_g, rhob_g, rhob_r)
740 CALL rs_grid_release(rden)
741 DEALLOCATE (rden)
742
743 DEALLOCATE (rhos)
744 DEALLOCATE (center, delta)
745 CALL timestop(handle)
746
747 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 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 NULLIFY (shell)
782 use_charge_array = PRESENT(charges)
783 IF (use_charge_array) use_charge_array = ASSOCIATED(charges)
784 IF (is_core .AND. is_shell) THEN
785 cpabort("Shell-model: cannot be core and shell simultaneously")
786 END IF
787
788 nbox = SIZE(rhos, 1)
789 r = part(p)%r
790 q = 1.0_dp
791 IF (.NOT. unit_charge) THEN
792 IF (is_core) THEN
793 CALL get_atomic_kind(atomic_kind=part(p)%atomic_kind, shell=shell)
794 q = shell%charge_core
795 ELSE IF (is_shell) THEN
796 CALL get_atomic_kind(atomic_kind=part(p)%atomic_kind, shell=shell)
797 q = shell%charge_shell
798 ELSE
799 CALL get_atomic_kind(atomic_kind=part(p)%atomic_kind, qeff=q)
800 END IF
801 IF (use_charge_array) q = charges(p)
802 END IF
803 CALL spme_get_patch(rhos, nbox, delta(:, p), q, green%p3m_coeff)
804
805 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 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 nbox = SIZE(rhos, 1)
833 r = part(p)%r
834 q = charges(p)
835 CALL spme_get_patch(rhos, nbox, delta(:, p), q, green%p3m_coeff)
836
837 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 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 IF (n > nmax) THEN
868 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 deltal(1, 0) = 1.0_dp
874 deltal(2, 0) = 1.0_dp
875 deltal(3, 0) = 1.0_dp
876 DO l = 1, n - 1
877 deltal(1, l) = deltal(1, l - 1)*delta(1)
878 deltal(2, l) = deltal(2, l - 1)*delta(2)
879 deltal(3, l) = deltal(3, l - 1)*delta(3)
880 END DO
881
882 w_assign = 0.0_dp
883 DO j = -(n - 1), n - 1, 2
884 DO l = 0, n - 1
885 w_assign(1, j) = w_assign(1, j) + coeff(j, l)*deltal(1, l)
886 w_assign(2, j) = w_assign(2, j) + coeff(j, l)*deltal(2, l)
887 w_assign(3, j) = w_assign(3, j) + coeff(j, l)*deltal(3, l)
888 END DO
889 END DO
890 DO i = 1, n
891 j = n + 1 - 2*i
892 f_assign(1, i) = w_assign(1, j)
893 f_assign(2, i) = w_assign(2, j)
894 f_assign(3, i) = w_assign(3, j)
895 END DO
896
897 DO i3 = 1, n
898 r3 = q*f_assign(3, i3)
899 DO i2 = 1, n
900 r2 = r3*f_assign(2, i2)
901 DO i1 = 1, n
902 rhos(i1, i2, i3) = r2*f_assign(1, i1)
903 END DO
904 END DO
905 END DO
906
907 END SUBROUTINE spme_get_patch
908
909END MODULE spme
910
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Holds information on atomic properties.
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public essmann1995
Handles all functions related to the CELL.
Definition cell_types.F:15
Definition dgs.F:13
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
subroutine, public ewald_pw_get(ewald_pw, pw_big_pool, pw_small_pool, rs_desc, poisson_env, dg)
get the ewald_pw environment to the correct program.
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Definition of mathematical constants and functions.
real(kind=dp), parameter, public fourpi
Interface to the message passing library MPI.
Define the data structure for the particle information.
Tools common both to PME and SPME.
Definition pme_tools.F:15
subroutine, public get_center(part, box, center, delta, npts, n)
...
Definition pme_tools.F:131
subroutine, public set_list(part, npart, center, p1, rs, ipart, core_center)
...
Definition pme_tools.F:46
This module defines the grid data type and some basic operations on it.
Definition pw_grids.F:36
subroutine, public get_pw_grid_info(pw_grid, id_nr, mode, vol, dvol, npts, ngpts, ngpts_cut, dr, cutoff, orthorhombic, gvectors, gsquare)
Access to information stored in the pw_grid_type.
Definition pw_grids.F:184
subroutine, public pw_derive(pw, n)
Calculate the derivative of a plane wave vector.
functions related to the poisson solver on regular grids
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
subroutine, public rs_grid_create(rs, desc)
...
subroutine, public transfer_pw2rs(rs, pw)
...
subroutine, public rs_grid_set_box(pw_grid, rs)
Set box matrix info for real space grid This is needed for variable cell simulations.
subroutine, public transfer_rs2pw(rs, pw)
...
subroutine, public rs_grid_release(rs_grid)
releases the given rs grid (see doc/ReferenceCounting.html)
subroutine, public rs_grid_zero(rs)
Initialize grid to zero.
Calculate the electrostatic energy by the Smooth Particle Ewald method.
Definition spme.F:14
subroutine, public spme_forces(ewald_env, ewald_pw, box, particle_set_a, charges_a, particle_set_b, charges_b, forces_b)
Calculate the forces on particles B for the electrostatic interaction betrween particles A and B.
Definition spme.F:608
subroutine, public spme_potential(ewald_env, ewald_pw, box, particle_set_a, charges_a, particle_set_b, potential)
Calculate the electrostatic potential from particles A (charge A) at positions of particles B.
Definition spme.F:467
subroutine, public spme_evaluate(ewald_env, ewald_pw, box, particle_set, fg_coulomb, vg_coulomb, pv_g, shell_particle_set, core_particle_set, fgshell_coulomb, fgcore_coulomb, use_virial, charges, atprop)
...
Definition spme.F:95
type for the atomic properties
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
contains all the informations needed by the fft based poisson solvers
environment for the poisson solver
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...