Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of Ewald contributions in DFTB
10 : !> \author JGH
11 : ! **************************************************************************************************
12 : MODULE ewald_methods_tb
13 : USE cell_types, ONLY: cell_type
14 : USE dgs, ONLY: dg_sum_patch,&
15 : dg_sum_patch_force_1d,&
16 : dg_sum_patch_force_3d
17 : USE ewald_environment_types, ONLY: ewald_env_get,&
18 : ewald_environment_type
19 : USE ewald_pw_types, ONLY: ewald_pw_get,&
20 : ewald_pw_type
21 : USE kinds, ONLY: dp
22 : USE mathconstants, ONLY: fourpi,&
23 : oorootpi
24 : USE message_passing, ONLY: mp_comm_type,&
25 : mp_para_env_type
26 : USE particle_types, ONLY: particle_type
27 : USE pme_tools, ONLY: get_center,&
28 : set_list
29 : USE pw_grid_types, ONLY: pw_grid_type
30 : USE pw_grids, ONLY: get_pw_grid_info
31 : USE pw_methods, ONLY: pw_integral_a2b,&
32 : pw_multiply_with,&
33 : pw_transfer
34 : USE pw_poisson_methods, ONLY: pw_poisson_rebuild,&
35 : pw_poisson_solve
36 : USE pw_poisson_types, ONLY: greens_fn_type,&
37 : pw_poisson_type
38 : USE pw_pool_types, ONLY: pw_pool_type
39 : USE pw_types, ONLY: pw_c1d_gs_type,&
40 : pw_r3d_rs_type
41 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
42 : neighbor_list_iterate,&
43 : neighbor_list_iterator_create,&
44 : neighbor_list_iterator_p_type,&
45 : neighbor_list_iterator_release,&
46 : neighbor_list_set_p_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 spme, ONLY: get_patch
56 : USE virial_methods, ONLY: virial_pair_force
57 : USE virial_types, ONLY: virial_type
58 : #include "./base/base_uses.f90"
59 :
60 : IMPLICIT NONE
61 :
62 : PRIVATE
63 :
64 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ewald_methods_tb'
65 :
66 : PUBLIC :: tb_spme_evaluate, tb_ewald_overlap, tb_spme_zforce
67 :
68 : CONTAINS
69 :
70 : ! **************************************************************************************************
71 : !> \brief ...
72 : !> \param ewald_env ...
73 : !> \param ewald_pw ...
74 : !> \param particle_set ...
75 : !> \param box ...
76 : !> \param gmcharge ...
77 : !> \param mcharge ...
78 : !> \param calculate_forces ...
79 : !> \param virial ...
80 : !> \param use_virial ...
81 : ! **************************************************************************************************
82 23522 : SUBROUTINE tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, &
83 23522 : gmcharge, mcharge, calculate_forces, virial, use_virial)
84 :
85 : TYPE(ewald_environment_type), POINTER :: ewald_env
86 : TYPE(ewald_pw_type), POINTER :: ewald_pw
87 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
88 : TYPE(cell_type), POINTER :: box
89 : REAL(KIND=dp), DIMENSION(:, :), INTENT(inout) :: gmcharge
90 : REAL(KIND=dp), DIMENSION(:), INTENT(in) :: mcharge
91 : LOGICAL, INTENT(in) :: calculate_forces
92 : TYPE(virial_type), POINTER :: virial
93 : LOGICAL, INTENT(in) :: use_virial
94 :
95 : CHARACTER(len=*), PARAMETER :: routineN = 'tb_spme_evaluate'
96 :
97 : INTEGER :: handle, i, ipart, j, n, npart, o_spline, &
98 : p1
99 23522 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: center
100 : INTEGER, DIMENSION(3) :: npts
101 : REAL(KIND=dp) :: alpha, dvols, fat(3), ffa, fint, vgc
102 23522 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: delta
103 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rhos
104 : REAL(KIND=dp), DIMENSION(3, 3) :: f_stress, h_stress
105 : TYPE(greens_fn_type), POINTER :: green
106 : TYPE(mp_comm_type) :: group
107 : TYPE(mp_para_env_type), POINTER :: para_env
108 94088 : TYPE(pw_c1d_gs_type), DIMENSION(3) :: dphi_g
109 : TYPE(pw_c1d_gs_type), POINTER :: phi_g, rhob_g
110 : TYPE(pw_grid_type), POINTER :: grid_spme
111 : TYPE(pw_poisson_type), POINTER :: poisson_env
112 : TYPE(pw_pool_type), POINTER :: pw_pool
113 : TYPE(pw_r3d_rs_type), POINTER :: rhob_r
114 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
115 917358 : TYPE(realspace_grid_type) :: rden, rpot
116 : TYPE(realspace_grid_type), ALLOCATABLE, &
117 23522 : DIMENSION(:) :: drpot
118 :
119 23522 : CALL timeset(routineN, handle)
120 : !-------------- INITIALISATION ---------------------
121 : CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group, &
122 23522 : para_env=para_env)
123 23522 : NULLIFY (green, poisson_env, pw_pool)
124 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, &
125 23522 : poisson_env=poisson_env)
126 23522 : CALL pw_poisson_rebuild(poisson_env)
127 23522 : green => poisson_env%green_fft
128 23522 : grid_spme => pw_pool%pw_grid
129 :
130 23522 : CALL get_pw_grid_info(grid_spme, dvol=dvols, npts=npts)
131 :
132 23522 : npart = SIZE(particle_set)
133 :
134 23522 : n = o_spline
135 117610 : ALLOCATE (rhos(n, n, n))
136 :
137 23522 : CALL rs_grid_create(rden, rs_desc)
138 23522 : CALL rs_grid_set_box(grid_spme, rs=rden)
139 23522 : CALL rs_grid_zero(rden)
140 :
141 117610 : ALLOCATE (center(3, npart), delta(3, npart))
142 23522 : CALL get_center(particle_set, box, center, delta, npts, n)
143 :
144 : !-------------- DENSITY CALCULATION ----------------
145 23522 : ipart = 0
146 145693 : DO
147 169215 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
148 169215 : IF (p1 == 0) EXIT
149 :
150 : ! calculate function on small boxes
151 : CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
152 145693 : is_shell=.FALSE., unit_charge=.TRUE.)
153 36721106 : rhos(:, :, :) = rhos(:, :, :)*mcharge(p1)
154 :
155 : ! add boxes to real space grid (big box)
156 169215 : CALL dg_sum_patch(rden, rhos, center(:, p1))
157 : END DO
158 :
159 23522 : NULLIFY (rhob_r)
160 23522 : ALLOCATE (rhob_r)
161 23522 : CALL pw_pool%create_pw(rhob_r)
162 :
163 23522 : CALL transfer_rs2pw(rden, rhob_r)
164 :
165 : ! transform density to G space and add charge function
166 23522 : NULLIFY (rhob_g)
167 23522 : ALLOCATE (rhob_g)
168 23522 : CALL pw_pool%create_pw(rhob_g)
169 23522 : CALL pw_transfer(rhob_r, rhob_g)
170 : ! update charge function
171 23522 : CALL pw_multiply_with(rhob_g, green%p3m_charge)
172 :
173 : !-------------- ELECTROSTATIC CALCULATION -----------
174 :
175 : ! allocate intermediate arrays
176 94088 : DO i = 1, 3
177 94088 : CALL pw_pool%create_pw(dphi_g(i))
178 : END DO
179 23522 : NULLIFY (phi_g)
180 23522 : ALLOCATE (phi_g)
181 23522 : CALL pw_pool%create_pw(phi_g)
182 23522 : IF (use_virial) THEN
183 224 : CALL pw_poisson_solve(poisson_env, rhob_g, vgc, phi_g, dphi_g, h_stress=h_stress)
184 : ELSE
185 23298 : CALL pw_poisson_solve(poisson_env, rhob_g, vgc, phi_g, dphi_g)
186 : END IF
187 :
188 23522 : CALL rs_grid_create(rpot, rs_desc)
189 23522 : CALL rs_grid_set_box(grid_spme, rs=rpot)
190 :
191 23522 : CALL pw_pool%give_back_pw(rhob_g)
192 23522 : DEALLOCATE (rhob_g)
193 :
194 23522 : CALL rs_grid_zero(rpot)
195 23522 : CALL pw_multiply_with(phi_g, green%p3m_charge)
196 23522 : CALL pw_transfer(phi_g, rhob_r)
197 23522 : CALL pw_pool%give_back_pw(phi_g)
198 23522 : DEALLOCATE (phi_g)
199 23522 : CALL transfer_pw2rs(rpot, rhob_r)
200 :
201 : !---------- END OF ELECTROSTATIC CALCULATION --------
202 :
203 : !------------- STRESS TENSOR CALCULATION ------------
204 :
205 23522 : IF (use_virial) THEN
206 896 : DO i = 1, 3
207 2240 : DO j = i, 3
208 1344 : f_stress(i, j) = pw_integral_a2b(dphi_g(i), dphi_g(j))
209 2016 : f_stress(j, i) = f_stress(i, j)
210 : END DO
211 : END DO
212 224 : ffa = (1.0_dp/fourpi)*(0.5_dp/alpha)**2
213 2912 : virial%pv_virial = virial%pv_virial - (ffa*f_stress - h_stress)/REAL(para_env%num_pe, dp)
214 : END IF
215 :
216 : !--------END OF STRESS TENSOR CALCULATION -----------
217 :
218 23522 : IF (calculate_forces) THEN
219 : ! move derivative of potential to real space grid and
220 : ! multiply by charge function in g-space
221 19044 : ALLOCATE (drpot(3))
222 3312 : DO i = 1, 3
223 2484 : CALL rs_grid_create(drpot(i), rs_desc)
224 2484 : CALL rs_grid_set_box(grid_spme, rs=drpot(i))
225 2484 : CALL pw_multiply_with(dphi_g(i), green%p3m_charge)
226 2484 : CALL pw_transfer(dphi_g(i), rhob_r)
227 2484 : CALL pw_pool%give_back_pw(dphi_g(i))
228 3312 : CALL transfer_pw2rs(drpot(i), rhob_r)
229 : END DO
230 : ELSE
231 90776 : DO i = 1, 3
232 90776 : CALL pw_pool%give_back_pw(dphi_g(i))
233 : END DO
234 : END IF
235 23522 : CALL pw_pool%give_back_pw(rhob_r)
236 23522 : DEALLOCATE (rhob_r)
237 :
238 : !----------------- FORCE CALCULATION ----------------
239 :
240 23522 : ipart = 0
241 : DO
242 :
243 169215 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
244 169215 : IF (p1 == 0) EXIT
245 :
246 : ! calculate function on small boxes
247 : CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
248 145693 : is_shell=.FALSE., unit_charge=.TRUE.)
249 :
250 145693 : CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fint)
251 145693 : gmcharge(p1, 1) = gmcharge(p1, 1) + fint*dvols
252 :
253 169215 : IF (calculate_forces) THEN
254 8585 : CALL dg_sum_patch_force_3d(drpot, rhos, center(:, p1), fat)
255 8585 : gmcharge(p1, 2) = gmcharge(p1, 2) - fat(1)*dvols
256 8585 : gmcharge(p1, 3) = gmcharge(p1, 3) - fat(2)*dvols
257 8585 : gmcharge(p1, 4) = gmcharge(p1, 4) - fat(3)*dvols
258 : END IF
259 :
260 : END DO
261 :
262 : !--------------END OF FORCE CALCULATION -------------
263 :
264 : !------------------CLEANING UP ----------------------
265 :
266 23522 : CALL rs_grid_release(rden)
267 23522 : CALL rs_grid_release(rpot)
268 23522 : IF (calculate_forces) THEN
269 3312 : DO i = 1, 3
270 3312 : CALL rs_grid_release(drpot(i))
271 : END DO
272 3312 : DEALLOCATE (drpot)
273 : END IF
274 23522 : DEALLOCATE (rhos)
275 23522 : DEALLOCATE (center, delta)
276 :
277 23522 : CALL timestop(handle)
278 :
279 70566 : END SUBROUTINE tb_spme_evaluate
280 :
281 : ! **************************************************************************************************
282 : !> \brief ...
283 : !> \param gmcharge ...
284 : !> \param mcharge ...
285 : !> \param alpha ...
286 : !> \param n_list ...
287 : !> \param virial ...
288 : !> \param use_virial ...
289 : ! **************************************************************************************************
290 21046 : SUBROUTINE tb_ewald_overlap(gmcharge, mcharge, alpha, n_list, virial, use_virial)
291 :
292 : REAL(KIND=dp), DIMENSION(:, :), INTENT(inout) :: gmcharge
293 : REAL(KIND=dp), DIMENSION(:), INTENT(in) :: mcharge
294 : REAL(KIND=dp), INTENT(in) :: alpha
295 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
296 : POINTER :: n_list
297 : TYPE(virial_type), POINTER :: virial
298 : LOGICAL, INTENT(IN) :: use_virial
299 :
300 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tb_ewald_overlap'
301 :
302 : INTEGER :: handle, i, iatom, jatom, nmat
303 : REAL(KIND=dp) :: dfr, dr, fr, pfr, rij(3)
304 : TYPE(neighbor_list_iterator_p_type), &
305 21046 : DIMENSION(:), POINTER :: nl_iterator
306 :
307 21046 : CALL timeset(routineN, handle)
308 :
309 21046 : nmat = SIZE(gmcharge, 2)
310 :
311 21046 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
312 133225700 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
313 133204654 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, r=rij)
314 :
315 532818616 : dr = SQRT(SUM(rij(:)**2))
316 133225700 : IF (dr > 1.e-10) THEN
317 133062494 : fr = erfc(alpha*dr)/dr
318 133062494 : gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge(jatom)*fr
319 133062494 : gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge(iatom)*fr
320 133062494 : IF (nmat > 1) THEN
321 11227759 : dfr = -2._dp*alpha*EXP(-alpha*alpha*dr*dr)*oorootpi/dr - fr/dr
322 11227759 : dfr = -dfr/dr
323 44911036 : DO i = 2, nmat
324 33683277 : gmcharge(iatom, i) = gmcharge(iatom, i) - rij(i - 1)*mcharge(jatom)*dfr
325 44911036 : gmcharge(jatom, i) = gmcharge(jatom, i) + rij(i - 1)*mcharge(iatom)*dfr
326 : END DO
327 : END IF
328 133062494 : IF (use_virial) THEN
329 : ! The potential update above adds both sides of the pair before the
330 : ! final 0.5 energy prefactor is applied. Self-image pairs therefore
331 : ! need the same virial prefactor as interatomic pairs.
332 7019681 : pfr = -dfr*mcharge(iatom)*mcharge(jatom)
333 7019681 : CALL virial_pair_force(virial%pv_virial, -pfr, rij, rij)
334 : END IF
335 : END IF
336 :
337 : END DO
338 21046 : CALL neighbor_list_iterator_release(nl_iterator)
339 :
340 21046 : CALL timestop(handle)
341 :
342 21046 : END SUBROUTINE tb_ewald_overlap
343 :
344 : ! **************************************************************************************************
345 : !> \brief ...
346 : !> \param ewald_env ...
347 : !> \param ewald_pw ...
348 : !> \param particle_set ...
349 : !> \param box ...
350 : !> \param gmcharge ...
351 : !> \param mcharge ...
352 : ! **************************************************************************************************
353 12 : SUBROUTINE tb_spme_zforce(ewald_env, ewald_pw, particle_set, box, gmcharge, mcharge)
354 :
355 : TYPE(ewald_environment_type), POINTER :: ewald_env
356 : TYPE(ewald_pw_type), POINTER :: ewald_pw
357 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
358 : TYPE(cell_type), POINTER :: box
359 : REAL(KIND=dp), DIMENSION(:, :), INTENT(inout) :: gmcharge
360 : REAL(KIND=dp), DIMENSION(:), INTENT(in) :: mcharge
361 :
362 : CHARACTER(len=*), PARAMETER :: routineN = 'tb_spme_zforce'
363 :
364 : INTEGER :: handle, i, ipart, n, npart, o_spline, p1
365 12 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: center
366 : INTEGER, DIMENSION(3) :: npts
367 : REAL(KIND=dp) :: alpha, dvols, fat(3), fint, vgc
368 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: delta
369 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rhos
370 : TYPE(greens_fn_type), POINTER :: green
371 : TYPE(mp_comm_type) :: group
372 : TYPE(mp_para_env_type), POINTER :: para_env
373 48 : TYPE(pw_c1d_gs_type), DIMENSION(3) :: dphi_g
374 : TYPE(pw_c1d_gs_type), POINTER :: phi_g, rhob_g
375 : TYPE(pw_grid_type), POINTER :: grid_spme
376 : TYPE(pw_poisson_type), POINTER :: poisson_env
377 : TYPE(pw_pool_type), POINTER :: pw_pool
378 : TYPE(pw_r3d_rs_type), POINTER :: rhob_r
379 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
380 468 : TYPE(realspace_grid_type) :: rden, rpot
381 432 : TYPE(realspace_grid_type), DIMENSION(3) :: drpot
382 :
383 12 : CALL timeset(routineN, handle)
384 : !-------------- INITIALISATION ---------------------
385 : CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group, &
386 12 : para_env=para_env)
387 12 : NULLIFY (green, poisson_env, pw_pool)
388 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, &
389 12 : poisson_env=poisson_env)
390 12 : CALL pw_poisson_rebuild(poisson_env)
391 12 : green => poisson_env%green_fft
392 12 : grid_spme => pw_pool%pw_grid
393 :
394 12 : CALL get_pw_grid_info(grid_spme, dvol=dvols, npts=npts)
395 :
396 12 : npart = SIZE(particle_set)
397 :
398 12 : n = o_spline
399 60 : ALLOCATE (rhos(n, n, n))
400 :
401 12 : CALL rs_grid_create(rden, rs_desc)
402 12 : CALL rs_grid_set_box(grid_spme, rs=rden)
403 12 : CALL rs_grid_zero(rden)
404 :
405 60 : ALLOCATE (center(3, npart), delta(3, npart))
406 12 : CALL get_center(particle_set, box, center, delta, npts, n)
407 :
408 : !-------------- DENSITY CALCULATION ----------------
409 12 : ipart = 0
410 206 : DO
411 218 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
412 218 : IF (p1 == 0) EXIT
413 :
414 : ! calculate function on small boxes
415 : CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
416 206 : is_shell=.FALSE., unit_charge=.TRUE.)
417 53354 : rhos(:, :, :) = rhos(:, :, :)*mcharge(p1)
418 :
419 : ! add boxes to real space grid (big box)
420 218 : CALL dg_sum_patch(rden, rhos, center(:, p1))
421 : END DO
422 :
423 12 : NULLIFY (rhob_r)
424 12 : ALLOCATE (rhob_r)
425 12 : CALL pw_pool%create_pw(rhob_r)
426 :
427 12 : CALL transfer_rs2pw(rden, rhob_r)
428 :
429 : ! transform density to G space and add charge function
430 12 : NULLIFY (rhob_g)
431 12 : ALLOCATE (rhob_g)
432 12 : CALL pw_pool%create_pw(rhob_g)
433 12 : CALL pw_transfer(rhob_r, rhob_g)
434 : ! update charge function
435 12 : CALL pw_multiply_with(rhob_g, green%p3m_charge)
436 :
437 : !-------------- ELECTROSTATIC CALCULATION -----------
438 :
439 : ! allocate intermediate arrays
440 48 : DO i = 1, 3
441 48 : CALL pw_pool%create_pw(dphi_g(i))
442 : END DO
443 12 : NULLIFY (phi_g)
444 12 : ALLOCATE (phi_g)
445 12 : CALL pw_pool%create_pw(phi_g)
446 12 : CALL pw_poisson_solve(poisson_env, rhob_g, vgc, phi_g, dphi_g)
447 :
448 12 : CALL rs_grid_create(rpot, rs_desc)
449 12 : CALL rs_grid_set_box(grid_spme, rs=rpot)
450 :
451 12 : CALL pw_pool%give_back_pw(rhob_g)
452 12 : DEALLOCATE (rhob_g)
453 :
454 12 : CALL rs_grid_zero(rpot)
455 12 : CALL pw_multiply_with(phi_g, green%p3m_charge)
456 12 : CALL pw_transfer(phi_g, rhob_r)
457 12 : CALL pw_pool%give_back_pw(phi_g)
458 12 : DEALLOCATE (phi_g)
459 12 : CALL transfer_pw2rs(rpot, rhob_r)
460 :
461 : !---------- END OF ELECTROSTATIC CALCULATION --------
462 :
463 : ! move derivative of potential to real space grid and
464 : ! multiply by charge function in g-space
465 48 : DO i = 1, 3
466 36 : CALL rs_grid_create(drpot(i), rs_desc)
467 36 : CALL rs_grid_set_box(grid_spme, rs=drpot(i))
468 36 : CALL pw_multiply_with(dphi_g(i), green%p3m_charge)
469 36 : CALL pw_transfer(dphi_g(i), rhob_r)
470 36 : CALL pw_pool%give_back_pw(dphi_g(i))
471 48 : CALL transfer_pw2rs(drpot(i), rhob_r)
472 : END DO
473 12 : CALL pw_pool%give_back_pw(rhob_r)
474 12 : DEALLOCATE (rhob_r)
475 :
476 : !----------------- FORCE CALCULATION ----------------
477 :
478 12 : ipart = 0
479 206 : DO
480 :
481 218 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
482 218 : IF (p1 == 0) EXIT
483 :
484 : ! calculate function on small boxes
485 : CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
486 206 : is_shell=.FALSE., unit_charge=.TRUE.)
487 :
488 206 : CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fint)
489 206 : gmcharge(p1, 1) = gmcharge(p1, 1) + fint*dvols
490 :
491 206 : CALL dg_sum_patch_force_3d(drpot, rhos, center(:, p1), fat)
492 206 : gmcharge(p1, 2) = gmcharge(p1, 2) - fat(1)*dvols
493 206 : gmcharge(p1, 3) = gmcharge(p1, 3) - fat(2)*dvols
494 206 : gmcharge(p1, 4) = gmcharge(p1, 4) - fat(3)*dvols
495 :
496 : END DO
497 :
498 : !--------------END OF FORCE CALCULATION -------------
499 :
500 : !------------------CLEANING UP ----------------------
501 :
502 12 : CALL rs_grid_release(rden)
503 12 : CALL rs_grid_release(rpot)
504 48 : DO i = 1, 3
505 48 : CALL rs_grid_release(drpot(i))
506 : END DO
507 12 : DEALLOCATE (rhos)
508 12 : DEALLOCATE (center, delta)
509 :
510 12 : CALL timestop(handle)
511 :
512 36 : END SUBROUTINE tb_spme_zforce
513 :
514 : END MODULE ewald_methods_tb
|
