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 atprop_types, ONLY: atprop_type
14 : USE cell_types, ONLY: cell_type
15 : USE dgs, ONLY: dg_sum_patch,&
16 : dg_sum_patch_force_1d,&
17 : dg_sum_patch_force_3d
18 : USE ewald_environment_types, ONLY: ewald_env_get,&
19 : ewald_environment_type
20 : USE ewald_pw_types, ONLY: ewald_pw_get,&
21 : ewald_pw_type
22 : USE kinds, ONLY: dp
23 : USE mathconstants, ONLY: fourpi,&
24 : oorootpi
25 : USE message_passing, ONLY: mp_comm_type,&
26 : mp_para_env_type
27 : USE particle_types, ONLY: particle_type
28 : USE pme_tools, ONLY: get_center,&
29 : set_list
30 : USE pw_grid_types, ONLY: pw_grid_type
31 : USE pw_grids, ONLY: get_pw_grid_info
32 : USE pw_methods, ONLY: pw_copy,&
33 : pw_derive,&
34 : pw_integral_a2b,&
35 : pw_multiply,&
36 : pw_multiply_with,&
37 : pw_transfer,&
38 : pw_zero
39 : USE pw_poisson_methods, ONLY: pw_poisson_rebuild,&
40 : pw_poisson_solve
41 : USE pw_poisson_types, ONLY: greens_fn_type,&
42 : pw_poisson_type
43 : USE pw_pool_types, ONLY: pw_pool_type
44 : USE pw_types, ONLY: pw_c1d_gs_type,&
45 : pw_r3d_rs_type
46 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
47 : neighbor_list_iterate,&
48 : neighbor_list_iterator_create,&
49 : neighbor_list_iterator_p_type,&
50 : neighbor_list_iterator_release,&
51 : neighbor_list_set_p_type
52 : USE realspace_grid_types, ONLY: realspace_grid_desc_type,&
53 : realspace_grid_type,&
54 : rs_grid_create,&
55 : rs_grid_release,&
56 : rs_grid_set_box,&
57 : rs_grid_zero,&
58 : transfer_pw2rs,&
59 : transfer_rs2pw
60 : USE spme, ONLY: get_patch
61 : USE virial_methods, ONLY: virial_pair_force
62 : USE virial_types, ONLY: virial_type
63 : #include "./base/base_uses.f90"
64 :
65 : IMPLICIT NONE
66 :
67 : PRIVATE
68 :
69 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ewald_methods_tb'
70 :
71 : PUBLIC :: tb_spme_evaluate, tb_ewald_overlap, tb_spme_zforce
72 :
73 : CONTAINS
74 :
75 : ! **************************************************************************************************
76 : !> \brief ...
77 : !> \param ewald_env ...
78 : !> \param ewald_pw ...
79 : !> \param particle_set ...
80 : !> \param box ...
81 : !> \param gmcharge ...
82 : !> \param mcharge ...
83 : !> \param calculate_forces ...
84 : !> \param virial ...
85 : !> \param use_virial ...
86 : !> \param atprop ...
87 : ! **************************************************************************************************
88 18408 : SUBROUTINE tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, &
89 18408 : gmcharge, mcharge, calculate_forces, virial, use_virial, atprop)
90 :
91 : TYPE(ewald_environment_type), POINTER :: ewald_env
92 : TYPE(ewald_pw_type), POINTER :: ewald_pw
93 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
94 : TYPE(cell_type), POINTER :: box
95 : REAL(KIND=dp), DIMENSION(:, :), INTENT(inout) :: gmcharge
96 : REAL(KIND=dp), DIMENSION(:), INTENT(in) :: mcharge
97 : LOGICAL, INTENT(in) :: calculate_forces
98 : TYPE(virial_type), POINTER :: virial
99 : LOGICAL, INTENT(in) :: use_virial
100 : TYPE(atprop_type), OPTIONAL, POINTER :: atprop
101 :
102 : CHARACTER(len=*), PARAMETER :: routineN = 'tb_spme_evaluate'
103 :
104 : INTEGER :: handle, i, ig, ipart, j, n, npart, &
105 : o_spline, p1
106 18408 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: center
107 : INTEGER, DIMENSION(3) :: nd, npts
108 : REAL(KIND=dp) :: alpha, dvols, fat(3), ffa, ffb, fint, vgc
109 18408 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: delta
110 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rhos
111 : REAL(KIND=dp), DIMENSION(3, 3) :: f_stress, h_stress
112 : TYPE(greens_fn_type), POINTER :: green
113 : TYPE(mp_comm_type) :: group
114 : TYPE(mp_para_env_type), POINTER :: para_env
115 73632 : TYPE(pw_c1d_gs_type), DIMENSION(3) :: dphi_g
116 : TYPE(pw_c1d_gs_type), POINTER :: phi_g, phib_g, rhob_g
117 : TYPE(pw_grid_type), POINTER :: grid_spme
118 : TYPE(pw_poisson_type), POINTER :: poisson_env
119 : TYPE(pw_pool_type), POINTER :: pw_pool
120 : TYPE(pw_r3d_rs_type), POINTER :: rhob_r
121 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
122 717912 : TYPE(realspace_grid_type) :: rden, rpot
123 : TYPE(realspace_grid_type), ALLOCATABLE, &
124 18408 : DIMENSION(:) :: drpot
125 :
126 18408 : CALL timeset(routineN, handle)
127 : !-------------- INITIALISATION ---------------------
128 : CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group, &
129 18408 : para_env=para_env)
130 18408 : NULLIFY (green, poisson_env, pw_pool)
131 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, &
132 18408 : poisson_env=poisson_env)
133 18408 : CALL pw_poisson_rebuild(poisson_env)
134 18408 : green => poisson_env%green_fft
135 18408 : grid_spme => pw_pool%pw_grid
136 :
137 18408 : CALL get_pw_grid_info(grid_spme, dvol=dvols, npts=npts)
138 :
139 18408 : npart = SIZE(particle_set)
140 :
141 18408 : n = o_spline
142 92040 : ALLOCATE (rhos(n, n, n))
143 :
144 18408 : CALL rs_grid_create(rden, rs_desc)
145 18408 : CALL rs_grid_set_box(grid_spme, rs=rden)
146 18408 : CALL rs_grid_zero(rden)
147 :
148 92040 : ALLOCATE (center(3, npart), delta(3, npart))
149 18408 : CALL get_center(particle_set, box, center, delta, npts, n)
150 :
151 : !-------------- DENSITY CALCULATION ----------------
152 18408 : ipart = 0
153 129125 : DO
154 147533 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
155 147533 : IF (p1 == 0) EXIT
156 :
157 : ! calculate function on small boxes
158 : CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
159 129125 : is_shell=.FALSE., unit_charge=.TRUE.)
160 32139225 : rhos(:, :, :) = rhos(:, :, :)*mcharge(p1)
161 :
162 : ! add boxes to real space grid (big box)
163 147533 : CALL dg_sum_patch(rden, rhos, center(:, p1))
164 : END DO
165 :
166 18408 : NULLIFY (rhob_r)
167 18408 : ALLOCATE (rhob_r)
168 18408 : CALL pw_pool%create_pw(rhob_r)
169 :
170 18408 : CALL transfer_rs2pw(rden, rhob_r)
171 :
172 : ! transform density to G space and add charge function
173 18408 : NULLIFY (rhob_g)
174 18408 : ALLOCATE (rhob_g)
175 18408 : CALL pw_pool%create_pw(rhob_g)
176 18408 : CALL pw_transfer(rhob_r, rhob_g)
177 : ! update charge function
178 18408 : CALL pw_multiply_with(rhob_g, green%p3m_charge)
179 :
180 : !-------------- ELECTROSTATIC CALCULATION -----------
181 :
182 : ! allocate intermediate arrays
183 73632 : DO i = 1, 3
184 73632 : CALL pw_pool%create_pw(dphi_g(i))
185 : END DO
186 18408 : NULLIFY (phi_g)
187 18408 : ALLOCATE (phi_g)
188 18408 : CALL pw_pool%create_pw(phi_g)
189 18408 : IF (use_virial) THEN
190 154 : CALL pw_poisson_solve(poisson_env, rhob_g, vgc, phi_g, dphi_g, h_stress=h_stress)
191 : ELSE
192 18254 : CALL pw_poisson_solve(poisson_env, rhob_g, vgc, phi_g, dphi_g)
193 : END IF
194 :
195 18408 : CALL rs_grid_create(rpot, rs_desc)
196 18408 : CALL rs_grid_set_box(grid_spme, rs=rpot)
197 :
198 : ! Atomic Stress
199 18408 : IF (PRESENT(atprop)) THEN
200 17490 : IF (ASSOCIATED(atprop)) THEN
201 17490 : IF (atprop%stress .AND. use_virial) THEN
202 :
203 52 : CALL rs_grid_zero(rpot)
204 52 : CALL pw_zero(rhob_g)
205 52 : CALL pw_multiply(rhob_g, phi_g, green%p3m_charge)
206 52 : CALL pw_transfer(rhob_g, rhob_r)
207 52 : CALL transfer_pw2rs(rpot, rhob_r)
208 52 : ipart = 0
209 2136 : DO
210 2188 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
211 2188 : IF (p1 == 0) EXIT
212 : ! calculate function on small boxes
213 : CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
214 2136 : is_shell=.FALSE., unit_charge=.TRUE.)
215 2136 : CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fint)
216 2136 : atprop%atstress(1, 1, p1) = atprop%atstress(1, 1, p1) + 0.5_dp*mcharge(p1)*fint*dvols
217 2136 : atprop%atstress(2, 2, p1) = atprop%atstress(2, 2, p1) + 0.5_dp*mcharge(p1)*fint*dvols
218 2136 : atprop%atstress(3, 3, p1) = atprop%atstress(3, 3, p1) + 0.5_dp*mcharge(p1)*fint*dvols
219 : END DO
220 :
221 52 : NULLIFY (phib_g)
222 52 : ALLOCATE (phib_g)
223 52 : CALL pw_pool%create_pw(phib_g)
224 52 : ffa = (0.5_dp/alpha)**2
225 52 : ffb = 1.0_dp/fourpi
226 208 : DO i = 1, 3
227 3741258 : DO ig = grid_spme%first_gne0, grid_spme%ngpts_cut_local
228 3741102 : phib_g%array(ig) = ffb*dphi_g(i)%array(ig)*(1.0_dp + ffa*grid_spme%gsq(ig))
229 3741258 : phib_g%array(ig) = phib_g%array(ig)*green%influence_fn%array(ig)
230 : END DO
231 156 : IF (grid_spme%have_g0) phib_g%array(1) = 0.0_dp
232 520 : DO j = 1, i
233 312 : nd = 0
234 312 : nd(j) = 1
235 312 : CALL pw_copy(phib_g, rhob_g)
236 312 : CALL pw_derive(rhob_g, nd)
237 312 : CALL pw_multiply_with(rhob_g, green%p3m_charge)
238 312 : CALL pw_transfer(rhob_g, rhob_r)
239 312 : CALL transfer_pw2rs(rpot, rhob_r)
240 :
241 312 : ipart = 0
242 156 : DO
243 13128 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
244 13128 : IF (p1 == 0) EXIT
245 : ! calculate function on small boxes
246 : CALL get_patch(particle_set, delta, green, p1, rhos, &
247 12816 : is_core=.FALSE., is_shell=.FALSE., unit_charge=.TRUE.)
248 : ! integrate box and potential
249 12816 : CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fint)
250 12816 : atprop%atstress(i, j, p1) = atprop%atstress(i, j, p1) + fint*dvols*mcharge(p1)
251 13128 : IF (i /= j) atprop%atstress(j, i, p1) = atprop%atstress(j, i, p1) + fint*dvols*mcharge(p1)
252 : END DO
253 :
254 : END DO
255 : END DO
256 52 : CALL pw_pool%give_back_pw(phib_g)
257 52 : DEALLOCATE (phib_g)
258 :
259 : END IF
260 : END IF
261 : END IF
262 :
263 18408 : CALL pw_pool%give_back_pw(rhob_g)
264 18408 : DEALLOCATE (rhob_g)
265 :
266 18408 : CALL rs_grid_zero(rpot)
267 18408 : CALL pw_multiply_with(phi_g, green%p3m_charge)
268 18408 : CALL pw_transfer(phi_g, rhob_r)
269 18408 : CALL pw_pool%give_back_pw(phi_g)
270 18408 : DEALLOCATE (phi_g)
271 18408 : CALL transfer_pw2rs(rpot, rhob_r)
272 :
273 : !---------- END OF ELECTROSTATIC CALCULATION --------
274 :
275 : !------------- STRESS TENSOR CALCULATION ------------
276 :
277 18408 : IF (use_virial) THEN
278 616 : DO i = 1, 3
279 1540 : DO j = i, 3
280 924 : f_stress(i, j) = pw_integral_a2b(dphi_g(i), dphi_g(j))
281 1386 : f_stress(j, i) = f_stress(i, j)
282 : END DO
283 : END DO
284 154 : ffa = (1.0_dp/fourpi)*(0.5_dp/alpha)**2
285 2002 : virial%pv_virial = virial%pv_virial - (ffa*f_stress - h_stress)/REAL(para_env%num_pe, dp)
286 : END IF
287 :
288 : !--------END OF STRESS TENSOR CALCULATION -----------
289 :
290 18408 : IF (calculate_forces) THEN
291 : ! move derivative of potential to real space grid and
292 : ! multiply by charge function in g-space
293 17434 : ALLOCATE (drpot(3))
294 3032 : DO i = 1, 3
295 2274 : CALL rs_grid_create(drpot(i), rs_desc)
296 2274 : CALL rs_grid_set_box(grid_spme, rs=drpot(i))
297 2274 : CALL pw_multiply_with(dphi_g(i), green%p3m_charge)
298 2274 : CALL pw_transfer(dphi_g(i), rhob_r)
299 2274 : CALL pw_pool%give_back_pw(dphi_g(i))
300 3032 : CALL transfer_pw2rs(drpot(i), rhob_r)
301 : END DO
302 : ELSE
303 70600 : DO i = 1, 3
304 70600 : CALL pw_pool%give_back_pw(dphi_g(i))
305 : END DO
306 : END IF
307 18408 : CALL pw_pool%give_back_pw(rhob_r)
308 18408 : DEALLOCATE (rhob_r)
309 :
310 : !----------------- FORCE CALCULATION ----------------
311 :
312 18408 : ipart = 0
313 : DO
314 :
315 147533 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
316 147533 : IF (p1 == 0) EXIT
317 :
318 : ! calculate function on small boxes
319 : CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
320 129125 : is_shell=.FALSE., unit_charge=.TRUE.)
321 :
322 129125 : CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fint)
323 129125 : gmcharge(p1, 1) = gmcharge(p1, 1) + fint*dvols
324 :
325 147533 : IF (calculate_forces) THEN
326 8330 : CALL dg_sum_patch_force_3d(drpot, rhos, center(:, p1), fat)
327 8330 : gmcharge(p1, 2) = gmcharge(p1, 2) - fat(1)*dvols
328 8330 : gmcharge(p1, 3) = gmcharge(p1, 3) - fat(2)*dvols
329 8330 : gmcharge(p1, 4) = gmcharge(p1, 4) - fat(3)*dvols
330 : END IF
331 :
332 : END DO
333 :
334 : !--------------END OF FORCE CALCULATION -------------
335 :
336 : !------------------CLEANING UP ----------------------
337 :
338 18408 : CALL rs_grid_release(rden)
339 18408 : CALL rs_grid_release(rpot)
340 18408 : IF (calculate_forces) THEN
341 3032 : DO i = 1, 3
342 3032 : CALL rs_grid_release(drpot(i))
343 : END DO
344 3032 : DEALLOCATE (drpot)
345 : END IF
346 18408 : DEALLOCATE (rhos)
347 18408 : DEALLOCATE (center, delta)
348 :
349 18408 : CALL timestop(handle)
350 :
351 55224 : END SUBROUTINE tb_spme_evaluate
352 :
353 : ! **************************************************************************************************
354 : !> \brief ...
355 : !> \param gmcharge ...
356 : !> \param mcharge ...
357 : !> \param alpha ...
358 : !> \param n_list ...
359 : !> \param virial ...
360 : !> \param use_virial ...
361 : !> \param atprop ...
362 : ! **************************************************************************************************
363 15932 : SUBROUTINE tb_ewald_overlap(gmcharge, mcharge, alpha, n_list, virial, use_virial, atprop)
364 :
365 : REAL(KIND=dp), DIMENSION(:, :), INTENT(inout) :: gmcharge
366 : REAL(KIND=dp), DIMENSION(:), INTENT(in) :: mcharge
367 : REAL(KIND=dp), INTENT(in) :: alpha
368 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
369 : POINTER :: n_list
370 : TYPE(virial_type), POINTER :: virial
371 : LOGICAL, INTENT(IN) :: use_virial
372 : TYPE(atprop_type), OPTIONAL, POINTER :: atprop
373 :
374 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tb_ewald_overlap'
375 :
376 : INTEGER :: handle, i, iatom, jatom, nmat
377 : REAL(KIND=dp) :: dfr, dr, fr, pfr, rij(3)
378 : TYPE(neighbor_list_iterator_p_type), &
379 15932 : DIMENSION(:), POINTER :: nl_iterator
380 :
381 15932 : CALL timeset(routineN, handle)
382 :
383 15932 : nmat = SIZE(gmcharge, 2)
384 :
385 15932 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
386 75718172 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
387 75702240 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, r=rij)
388 :
389 302808960 : dr = SQRT(SUM(rij(:)**2))
390 75718172 : IF (dr > 1.e-10) THEN
391 75576648 : fr = erfc(alpha*dr)/dr
392 75576648 : gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge(jatom)*fr
393 75576648 : gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge(iatom)*fr
394 75576648 : IF (nmat > 1) THEN
395 9359599 : dfr = -2._dp*alpha*EXP(-alpha*alpha*dr*dr)*oorootpi/dr - fr/dr
396 9359599 : dfr = -dfr/dr
397 37438396 : DO i = 2, nmat
398 28078797 : gmcharge(iatom, i) = gmcharge(iatom, i) - rij(i - 1)*mcharge(jatom)*dfr
399 37438396 : gmcharge(jatom, i) = gmcharge(jatom, i) + rij(i - 1)*mcharge(iatom)*dfr
400 : END DO
401 : END IF
402 75576648 : IF (use_virial) THEN
403 6231309 : IF (iatom == jatom) THEN
404 132432 : pfr = -0.5_dp*dfr*mcharge(iatom)*mcharge(jatom)
405 : ELSE
406 6098877 : pfr = -dfr*mcharge(iatom)*mcharge(jatom)
407 : END IF
408 6231309 : CALL virial_pair_force(virial%pv_virial, -pfr, rij, rij)
409 6231309 : IF (PRESENT(atprop)) THEN
410 6231309 : IF (ASSOCIATED(atprop)) THEN
411 6231309 : IF (atprop%stress) THEN
412 2620724 : CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp*pfr, rij, rij)
413 2620724 : CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp*pfr, rij, rij)
414 : END IF
415 : END IF
416 : END IF
417 : END IF
418 : END IF
419 :
420 : END DO
421 15932 : CALL neighbor_list_iterator_release(nl_iterator)
422 :
423 15932 : CALL timestop(handle)
424 :
425 15932 : END SUBROUTINE tb_ewald_overlap
426 :
427 : ! **************************************************************************************************
428 : !> \brief ...
429 : !> \param ewald_env ...
430 : !> \param ewald_pw ...
431 : !> \param particle_set ...
432 : !> \param box ...
433 : !> \param gmcharge ...
434 : !> \param mcharge ...
435 : ! **************************************************************************************************
436 12 : SUBROUTINE tb_spme_zforce(ewald_env, ewald_pw, particle_set, box, gmcharge, mcharge)
437 :
438 : TYPE(ewald_environment_type), POINTER :: ewald_env
439 : TYPE(ewald_pw_type), POINTER :: ewald_pw
440 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
441 : TYPE(cell_type), POINTER :: box
442 : REAL(KIND=dp), DIMENSION(:, :), INTENT(inout) :: gmcharge
443 : REAL(KIND=dp), DIMENSION(:), INTENT(in) :: mcharge
444 :
445 : CHARACTER(len=*), PARAMETER :: routineN = 'tb_spme_zforce'
446 :
447 : INTEGER :: handle, i, ipart, n, npart, o_spline, p1
448 12 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: center
449 : INTEGER, DIMENSION(3) :: npts
450 : REAL(KIND=dp) :: alpha, dvols, fat(3), fint, vgc
451 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: delta
452 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rhos
453 : TYPE(greens_fn_type), POINTER :: green
454 : TYPE(mp_comm_type) :: group
455 : TYPE(mp_para_env_type), POINTER :: para_env
456 48 : TYPE(pw_c1d_gs_type), DIMENSION(3) :: dphi_g
457 : TYPE(pw_c1d_gs_type), POINTER :: phi_g, rhob_g
458 : TYPE(pw_grid_type), POINTER :: grid_spme
459 : TYPE(pw_poisson_type), POINTER :: poisson_env
460 : TYPE(pw_pool_type), POINTER :: pw_pool
461 : TYPE(pw_r3d_rs_type), POINTER :: rhob_r
462 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
463 468 : TYPE(realspace_grid_type) :: rden, rpot
464 432 : TYPE(realspace_grid_type), DIMENSION(3) :: drpot
465 :
466 12 : CALL timeset(routineN, handle)
467 : !-------------- INITIALISATION ---------------------
468 : CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group, &
469 12 : para_env=para_env)
470 12 : NULLIFY (green, poisson_env, pw_pool)
471 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, &
472 12 : poisson_env=poisson_env)
473 12 : CALL pw_poisson_rebuild(poisson_env)
474 12 : green => poisson_env%green_fft
475 12 : grid_spme => pw_pool%pw_grid
476 :
477 12 : CALL get_pw_grid_info(grid_spme, dvol=dvols, npts=npts)
478 :
479 12 : npart = SIZE(particle_set)
480 :
481 12 : n = o_spline
482 60 : ALLOCATE (rhos(n, n, n))
483 :
484 12 : CALL rs_grid_create(rden, rs_desc)
485 12 : CALL rs_grid_set_box(grid_spme, rs=rden)
486 12 : CALL rs_grid_zero(rden)
487 :
488 60 : ALLOCATE (center(3, npart), delta(3, npart))
489 12 : CALL get_center(particle_set, box, center, delta, npts, n)
490 :
491 : !-------------- DENSITY CALCULATION ----------------
492 12 : ipart = 0
493 206 : DO
494 218 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
495 218 : IF (p1 == 0) EXIT
496 :
497 : ! calculate function on small boxes
498 : CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
499 206 : is_shell=.FALSE., unit_charge=.TRUE.)
500 53354 : rhos(:, :, :) = rhos(:, :, :)*mcharge(p1)
501 :
502 : ! add boxes to real space grid (big box)
503 218 : CALL dg_sum_patch(rden, rhos, center(:, p1))
504 : END DO
505 :
506 12 : NULLIFY (rhob_r)
507 12 : ALLOCATE (rhob_r)
508 12 : CALL pw_pool%create_pw(rhob_r)
509 :
510 12 : CALL transfer_rs2pw(rden, rhob_r)
511 :
512 : ! transform density to G space and add charge function
513 12 : NULLIFY (rhob_g)
514 12 : ALLOCATE (rhob_g)
515 12 : CALL pw_pool%create_pw(rhob_g)
516 12 : CALL pw_transfer(rhob_r, rhob_g)
517 : ! update charge function
518 12 : CALL pw_multiply_with(rhob_g, green%p3m_charge)
519 :
520 : !-------------- ELECTROSTATIC CALCULATION -----------
521 :
522 : ! allocate intermediate arrays
523 48 : DO i = 1, 3
524 48 : CALL pw_pool%create_pw(dphi_g(i))
525 : END DO
526 12 : NULLIFY (phi_g)
527 12 : ALLOCATE (phi_g)
528 12 : CALL pw_pool%create_pw(phi_g)
529 12 : CALL pw_poisson_solve(poisson_env, rhob_g, vgc, phi_g, dphi_g)
530 :
531 12 : CALL rs_grid_create(rpot, rs_desc)
532 12 : CALL rs_grid_set_box(grid_spme, rs=rpot)
533 :
534 12 : CALL pw_pool%give_back_pw(rhob_g)
535 12 : DEALLOCATE (rhob_g)
536 :
537 12 : CALL rs_grid_zero(rpot)
538 12 : CALL pw_multiply_with(phi_g, green%p3m_charge)
539 12 : CALL pw_transfer(phi_g, rhob_r)
540 12 : CALL pw_pool%give_back_pw(phi_g)
541 12 : DEALLOCATE (phi_g)
542 12 : CALL transfer_pw2rs(rpot, rhob_r)
543 :
544 : !---------- END OF ELECTROSTATIC CALCULATION --------
545 :
546 : ! move derivative of potential to real space grid and
547 : ! multiply by charge function in g-space
548 48 : DO i = 1, 3
549 36 : CALL rs_grid_create(drpot(i), rs_desc)
550 36 : CALL rs_grid_set_box(grid_spme, rs=drpot(i))
551 36 : CALL pw_multiply_with(dphi_g(i), green%p3m_charge)
552 36 : CALL pw_transfer(dphi_g(i), rhob_r)
553 36 : CALL pw_pool%give_back_pw(dphi_g(i))
554 48 : CALL transfer_pw2rs(drpot(i), rhob_r)
555 : END DO
556 12 : CALL pw_pool%give_back_pw(rhob_r)
557 12 : DEALLOCATE (rhob_r)
558 :
559 : !----------------- FORCE CALCULATION ----------------
560 :
561 12 : ipart = 0
562 206 : DO
563 :
564 218 : CALL set_list(particle_set, npart, center, p1, rden, ipart)
565 218 : IF (p1 == 0) EXIT
566 :
567 : ! calculate function on small boxes
568 : CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
569 206 : is_shell=.FALSE., unit_charge=.TRUE.)
570 :
571 206 : CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fint)
572 206 : gmcharge(p1, 1) = gmcharge(p1, 1) + fint*dvols
573 :
574 206 : CALL dg_sum_patch_force_3d(drpot, rhos, center(:, p1), fat)
575 206 : gmcharge(p1, 2) = gmcharge(p1, 2) - fat(1)*dvols
576 206 : gmcharge(p1, 3) = gmcharge(p1, 3) - fat(2)*dvols
577 206 : gmcharge(p1, 4) = gmcharge(p1, 4) - fat(3)*dvols
578 :
579 : END DO
580 :
581 : !--------------END OF FORCE CALCULATION -------------
582 :
583 : !------------------CLEANING UP ----------------------
584 :
585 12 : CALL rs_grid_release(rden)
586 12 : CALL rs_grid_release(rpot)
587 48 : DO i = 1, 3
588 48 : CALL rs_grid_release(drpot(i))
589 : END DO
590 12 : DEALLOCATE (rhos)
591 12 : DEALLOCATE (center, delta)
592 :
593 12 : CALL timestop(handle)
594 :
595 36 : END SUBROUTINE tb_spme_zforce
596 :
597 : END MODULE ewald_methods_tb
598 :
|
