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 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 18432 : SUBROUTINE tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, &
83 18432 : 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 18432 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: center
100 : INTEGER, DIMENSION(3) :: npts
101 : REAL(KIND=dp) :: alpha, dvols, fat(3), ffa, fint, vgc
102 18432 : 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 73728 : 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 718848 : TYPE(realspace_grid_type) :: rden, rpot
116 : TYPE(realspace_grid_type), ALLOCATABLE, &
117 18432 : DIMENSION(:) :: drpot
118 :
119 18432 : CALL timeset(routineN, handle)
120 : !-------------- INITIALISATION ---------------------
121 : CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group, &
122 18432 : para_env=para_env)
123 18432 : NULLIFY (green, poisson_env, pw_pool)
124 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, &
125 18432 : poisson_env=poisson_env)
126 18432 : CALL pw_poisson_rebuild(poisson_env)
127 18432 : green => poisson_env%green_fft
128 18432 : grid_spme => pw_pool%pw_grid
129 :
130 18432 : CALL get_pw_grid_info(grid_spme, dvol=dvols, npts=npts)
131 :
132 18432 : npart = SIZE(particle_set)
133 :
134 18432 : n = o_spline
135 92160 : ALLOCATE (rhos(n, n, n))
136 :
137 18432 : CALL rs_grid_create(rden, rs_desc)
138 18432 : CALL rs_grid_set_box(grid_spme, rs=rden)
139 18432 : CALL rs_grid_zero(rden)
140 :
141 92160 : ALLOCATE (center(3, npart), delta(3, npart))
142 18432 : CALL get_center(particle_set, box, center, delta, npts, n)
143 :
144 : !-------------- DENSITY CALCULATION ----------------
145 18432 : ipart = 0
146 128987 : DO
147 147419 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
148 147419 : 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 128987 : is_shell=.FALSE., unit_charge=.TRUE.)
153 32093595 : rhos(:, :, :) = rhos(:, :, :)*mcharge(p1)
154 :
155 : ! add boxes to real space grid (big box)
156 147419 : CALL dg_sum_patch(rden, rhos, center(:, p1))
157 : END DO
158 :
159 18432 : NULLIFY (rhob_r)
160 18432 : ALLOCATE (rhob_r)
161 18432 : CALL pw_pool%create_pw(rhob_r)
162 :
163 18432 : CALL transfer_rs2pw(rden, rhob_r)
164 :
165 : ! transform density to G space and add charge function
166 18432 : NULLIFY (rhob_g)
167 18432 : ALLOCATE (rhob_g)
168 18432 : CALL pw_pool%create_pw(rhob_g)
169 18432 : CALL pw_transfer(rhob_r, rhob_g)
170 : ! update charge function
171 18432 : CALL pw_multiply_with(rhob_g, green%p3m_charge)
172 :
173 : !-------------- ELECTROSTATIC CALCULATION -----------
174 :
175 : ! allocate intermediate arrays
176 73728 : DO i = 1, 3
177 73728 : CALL pw_pool%create_pw(dphi_g(i))
178 : END DO
179 18432 : NULLIFY (phi_g)
180 18432 : ALLOCATE (phi_g)
181 18432 : CALL pw_pool%create_pw(phi_g)
182 18432 : IF (use_virial) THEN
183 154 : CALL pw_poisson_solve(poisson_env, rhob_g, vgc, phi_g, dphi_g, h_stress=h_stress)
184 : ELSE
185 18278 : CALL pw_poisson_solve(poisson_env, rhob_g, vgc, phi_g, dphi_g)
186 : END IF
187 :
188 18432 : CALL rs_grid_create(rpot, rs_desc)
189 18432 : CALL rs_grid_set_box(grid_spme, rs=rpot)
190 :
191 18432 : CALL pw_pool%give_back_pw(rhob_g)
192 18432 : DEALLOCATE (rhob_g)
193 :
194 18432 : CALL rs_grid_zero(rpot)
195 18432 : CALL pw_multiply_with(phi_g, green%p3m_charge)
196 18432 : CALL pw_transfer(phi_g, rhob_r)
197 18432 : CALL pw_pool%give_back_pw(phi_g)
198 18432 : DEALLOCATE (phi_g)
199 18432 : CALL transfer_pw2rs(rpot, rhob_r)
200 :
201 : !---------- END OF ELECTROSTATIC CALCULATION --------
202 :
203 : !------------- STRESS TENSOR CALCULATION ------------
204 :
205 18432 : IF (use_virial) THEN
206 616 : DO i = 1, 3
207 1540 : DO j = i, 3
208 924 : f_stress(i, j) = pw_integral_a2b(dphi_g(i), dphi_g(j))
209 1386 : f_stress(j, i) = f_stress(i, j)
210 : END DO
211 : END DO
212 154 : ffa = (1.0_dp/fourpi)*(0.5_dp/alpha)**2
213 2002 : 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 18432 : IF (calculate_forces) THEN
219 : ! move derivative of potential to real space grid and
220 : ! multiply by charge function in g-space
221 17434 : ALLOCATE (drpot(3))
222 3032 : DO i = 1, 3
223 2274 : CALL rs_grid_create(drpot(i), rs_desc)
224 2274 : CALL rs_grid_set_box(grid_spme, rs=drpot(i))
225 2274 : CALL pw_multiply_with(dphi_g(i), green%p3m_charge)
226 2274 : CALL pw_transfer(dphi_g(i), rhob_r)
227 2274 : CALL pw_pool%give_back_pw(dphi_g(i))
228 3032 : CALL transfer_pw2rs(drpot(i), rhob_r)
229 : END DO
230 : ELSE
231 70696 : DO i = 1, 3
232 70696 : CALL pw_pool%give_back_pw(dphi_g(i))
233 : END DO
234 : END IF
235 18432 : CALL pw_pool%give_back_pw(rhob_r)
236 18432 : DEALLOCATE (rhob_r)
237 :
238 : !----------------- FORCE CALCULATION ----------------
239 :
240 18432 : ipart = 0
241 : DO
242 :
243 147419 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
244 147419 : 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 128987 : is_shell=.FALSE., unit_charge=.TRUE.)
249 :
250 128987 : CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fint)
251 128987 : gmcharge(p1, 1) = gmcharge(p1, 1) + fint*dvols
252 :
253 147419 : IF (calculate_forces) THEN
254 8330 : CALL dg_sum_patch_force_3d(drpot, rhos, center(:, p1), fat)
255 8330 : gmcharge(p1, 2) = gmcharge(p1, 2) - fat(1)*dvols
256 8330 : gmcharge(p1, 3) = gmcharge(p1, 3) - fat(2)*dvols
257 8330 : 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 18432 : CALL rs_grid_release(rden)
267 18432 : CALL rs_grid_release(rpot)
268 18432 : IF (calculate_forces) THEN
269 3032 : DO i = 1, 3
270 3032 : CALL rs_grid_release(drpot(i))
271 : END DO
272 3032 : DEALLOCATE (drpot)
273 : END IF
274 18432 : DEALLOCATE (rhos)
275 18432 : DEALLOCATE (center, delta)
276 :
277 18432 : CALL timestop(handle)
278 :
279 55296 : 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 15956 : 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 15956 : DIMENSION(:), POINTER :: nl_iterator
306 :
307 15956 : CALL timeset(routineN, handle)
308 :
309 15956 : nmat = SIZE(gmcharge, 2)
310 :
311 15956 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
312 75290922 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
313 75274966 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, r=rij)
314 :
315 301099864 : dr = SQRT(SUM(rij(:)**2))
316 75290922 : IF (dr > 1.e-10) THEN
317 75149512 : fr = erfc(alpha*dr)/dr
318 75149512 : gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge(jatom)*fr
319 75149512 : gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge(iatom)*fr
320 75149512 : IF (nmat > 1) THEN
321 9359599 : dfr = -2._dp*alpha*EXP(-alpha*alpha*dr*dr)*oorootpi/dr - fr/dr
322 9359599 : dfr = -dfr/dr
323 37438396 : DO i = 2, nmat
324 28078797 : gmcharge(iatom, i) = gmcharge(iatom, i) - rij(i - 1)*mcharge(jatom)*dfr
325 37438396 : gmcharge(jatom, i) = gmcharge(jatom, i) + rij(i - 1)*mcharge(iatom)*dfr
326 : END DO
327 : END IF
328 75149512 : IF (use_virial) THEN
329 6231309 : IF (iatom == jatom) THEN
330 132432 : pfr = -0.5_dp*dfr*mcharge(iatom)*mcharge(jatom)
331 : ELSE
332 6098877 : pfr = -dfr*mcharge(iatom)*mcharge(jatom)
333 : END IF
334 6231309 : CALL virial_pair_force(virial%pv_virial, -pfr, rij, rij)
335 : END IF
336 : END IF
337 :
338 : END DO
339 15956 : CALL neighbor_list_iterator_release(nl_iterator)
340 :
341 15956 : CALL timestop(handle)
342 :
343 15956 : END SUBROUTINE tb_ewald_overlap
344 :
345 : ! **************************************************************************************************
346 : !> \brief ...
347 : !> \param ewald_env ...
348 : !> \param ewald_pw ...
349 : !> \param particle_set ...
350 : !> \param box ...
351 : !> \param gmcharge ...
352 : !> \param mcharge ...
353 : ! **************************************************************************************************
354 12 : SUBROUTINE tb_spme_zforce(ewald_env, ewald_pw, particle_set, box, gmcharge, mcharge)
355 :
356 : TYPE(ewald_environment_type), POINTER :: ewald_env
357 : TYPE(ewald_pw_type), POINTER :: ewald_pw
358 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
359 : TYPE(cell_type), POINTER :: box
360 : REAL(KIND=dp), DIMENSION(:, :), INTENT(inout) :: gmcharge
361 : REAL(KIND=dp), DIMENSION(:), INTENT(in) :: mcharge
362 :
363 : CHARACTER(len=*), PARAMETER :: routineN = 'tb_spme_zforce'
364 :
365 : INTEGER :: handle, i, ipart, n, npart, o_spline, p1
366 12 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: center
367 : INTEGER, DIMENSION(3) :: npts
368 : REAL(KIND=dp) :: alpha, dvols, fat(3), fint, vgc
369 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: delta
370 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rhos
371 : TYPE(greens_fn_type), POINTER :: green
372 : TYPE(mp_comm_type) :: group
373 : TYPE(mp_para_env_type), POINTER :: para_env
374 48 : TYPE(pw_c1d_gs_type), DIMENSION(3) :: dphi_g
375 : TYPE(pw_c1d_gs_type), POINTER :: phi_g, rhob_g
376 : TYPE(pw_grid_type), POINTER :: grid_spme
377 : TYPE(pw_poisson_type), POINTER :: poisson_env
378 : TYPE(pw_pool_type), POINTER :: pw_pool
379 : TYPE(pw_r3d_rs_type), POINTER :: rhob_r
380 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
381 468 : TYPE(realspace_grid_type) :: rden, rpot
382 432 : TYPE(realspace_grid_type), DIMENSION(3) :: drpot
383 :
384 12 : CALL timeset(routineN, handle)
385 : !-------------- INITIALISATION ---------------------
386 : CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group, &
387 12 : para_env=para_env)
388 12 : NULLIFY (green, poisson_env, pw_pool)
389 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, &
390 12 : poisson_env=poisson_env)
391 12 : CALL pw_poisson_rebuild(poisson_env)
392 12 : green => poisson_env%green_fft
393 12 : grid_spme => pw_pool%pw_grid
394 :
395 12 : CALL get_pw_grid_info(grid_spme, dvol=dvols, npts=npts)
396 :
397 12 : npart = SIZE(particle_set)
398 :
399 12 : n = o_spline
400 60 : ALLOCATE (rhos(n, n, n))
401 :
402 12 : CALL rs_grid_create(rden, rs_desc)
403 12 : CALL rs_grid_set_box(grid_spme, rs=rden)
404 12 : CALL rs_grid_zero(rden)
405 :
406 60 : ALLOCATE (center(3, npart), delta(3, npart))
407 12 : CALL get_center(particle_set, box, center, delta, npts, n)
408 :
409 : !-------------- DENSITY CALCULATION ----------------
410 12 : ipart = 0
411 206 : DO
412 218 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
413 218 : IF (p1 == 0) EXIT
414 :
415 : ! calculate function on small boxes
416 : CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
417 206 : is_shell=.FALSE., unit_charge=.TRUE.)
418 53354 : rhos(:, :, :) = rhos(:, :, :)*mcharge(p1)
419 :
420 : ! add boxes to real space grid (big box)
421 218 : CALL dg_sum_patch(rden, rhos, center(:, p1))
422 : END DO
423 :
424 12 : NULLIFY (rhob_r)
425 12 : ALLOCATE (rhob_r)
426 12 : CALL pw_pool%create_pw(rhob_r)
427 :
428 12 : CALL transfer_rs2pw(rden, rhob_r)
429 :
430 : ! transform density to G space and add charge function
431 12 : NULLIFY (rhob_g)
432 12 : ALLOCATE (rhob_g)
433 12 : CALL pw_pool%create_pw(rhob_g)
434 12 : CALL pw_transfer(rhob_r, rhob_g)
435 : ! update charge function
436 12 : CALL pw_multiply_with(rhob_g, green%p3m_charge)
437 :
438 : !-------------- ELECTROSTATIC CALCULATION -----------
439 :
440 : ! allocate intermediate arrays
441 48 : DO i = 1, 3
442 48 : CALL pw_pool%create_pw(dphi_g(i))
443 : END DO
444 12 : NULLIFY (phi_g)
445 12 : ALLOCATE (phi_g)
446 12 : CALL pw_pool%create_pw(phi_g)
447 12 : CALL pw_poisson_solve(poisson_env, rhob_g, vgc, phi_g, dphi_g)
448 :
449 12 : CALL rs_grid_create(rpot, rs_desc)
450 12 : CALL rs_grid_set_box(grid_spme, rs=rpot)
451 :
452 12 : CALL pw_pool%give_back_pw(rhob_g)
453 12 : DEALLOCATE (rhob_g)
454 :
455 12 : CALL rs_grid_zero(rpot)
456 12 : CALL pw_multiply_with(phi_g, green%p3m_charge)
457 12 : CALL pw_transfer(phi_g, rhob_r)
458 12 : CALL pw_pool%give_back_pw(phi_g)
459 12 : DEALLOCATE (phi_g)
460 12 : CALL transfer_pw2rs(rpot, rhob_r)
461 :
462 : !---------- END OF ELECTROSTATIC CALCULATION --------
463 :
464 : ! move derivative of potential to real space grid and
465 : ! multiply by charge function in g-space
466 48 : DO i = 1, 3
467 36 : CALL rs_grid_create(drpot(i), rs_desc)
468 36 : CALL rs_grid_set_box(grid_spme, rs=drpot(i))
469 36 : CALL pw_multiply_with(dphi_g(i), green%p3m_charge)
470 36 : CALL pw_transfer(dphi_g(i), rhob_r)
471 36 : CALL pw_pool%give_back_pw(dphi_g(i))
472 48 : CALL transfer_pw2rs(drpot(i), rhob_r)
473 : END DO
474 12 : CALL pw_pool%give_back_pw(rhob_r)
475 12 : DEALLOCATE (rhob_r)
476 :
477 : !----------------- FORCE CALCULATION ----------------
478 :
479 12 : ipart = 0
480 206 : DO
481 :
482 218 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
483 218 : IF (p1 == 0) EXIT
484 :
485 : ! calculate function on small boxes
486 : CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
487 206 : is_shell=.FALSE., unit_charge=.TRUE.)
488 :
489 206 : CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fint)
490 206 : gmcharge(p1, 1) = gmcharge(p1, 1) + fint*dvols
491 :
492 206 : CALL dg_sum_patch_force_3d(drpot, rhos, center(:, p1), fat)
493 206 : gmcharge(p1, 2) = gmcharge(p1, 2) - fat(1)*dvols
494 206 : gmcharge(p1, 3) = gmcharge(p1, 3) - fat(2)*dvols
495 206 : gmcharge(p1, 4) = gmcharge(p1, 4) - fat(3)*dvols
496 :
497 : END DO
498 :
499 : !--------------END OF FORCE CALCULATION -------------
500 :
501 : !------------------CLEANING UP ----------------------
502 :
503 12 : CALL rs_grid_release(rden)
504 12 : CALL rs_grid_release(rpot)
505 48 : DO i = 1, 3
506 48 : CALL rs_grid_release(drpot(i))
507 : END DO
508 12 : DEALLOCATE (rhos)
509 12 : DEALLOCATE (center, delta)
510 :
511 12 : CALL timestop(handle)
512 :
513 36 : END SUBROUTINE tb_spme_zforce
514 :
515 : END MODULE ewald_methods_tb
516 :
|
