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 : !> \par History
10 : !> JGH (21-Mar-2001) : Complete rewrite
11 : !> \author CJM and APSI
12 : ! **************************************************************************************************
13 : MODULE pme
14 :
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind
17 : USE atprop_types, ONLY: atprop_type
18 : USE bibliography, ONLY: cite_reference,&
19 : darden1993
20 : USE cell_types, ONLY: cell_type
21 : USE dg_rho0_types, ONLY: dg_rho0_type
22 : USE dg_types, ONLY: dg_get,&
23 : dg_type
24 : USE dgs, ONLY: dg_get_patch,&
25 : dg_get_strucfac,&
26 : dg_sum_patch,&
27 : dg_sum_patch_force_1d,&
28 : dg_sum_patch_force_3d
29 : USE ewald_environment_types, ONLY: ewald_env_get,&
30 : ewald_environment_type
31 : USE ewald_pw_types, ONLY: ewald_pw_get,&
32 : ewald_pw_type
33 : USE kinds, ONLY: dp
34 : USE mathconstants, ONLY: fourpi
35 : USE message_passing, ONLY: mp_comm_type
36 : USE particle_types, ONLY: particle_type
37 : USE pme_tools, ONLY: get_center,&
38 : set_list
39 : USE pw_grid_types, ONLY: pw_grid_type
40 : USE pw_methods, ONLY: pw_copy,&
41 : pw_derive,&
42 : pw_integral_a2b,&
43 : pw_transfer
44 : USE pw_poisson_methods, ONLY: pw_poisson_solve
45 : USE pw_poisson_types, ONLY: pw_poisson_type
46 : USE pw_pool_types, ONLY: pw_pool_type
47 : USE pw_types, ONLY: pw_c1d_gs_type,&
48 : pw_r3d_rs_type
49 : USE realspace_grid_types, ONLY: realspace_grid_desc_type,&
50 : realspace_grid_type,&
51 : rs_grid_create,&
52 : rs_grid_release,&
53 : rs_grid_set_box,&
54 : rs_grid_zero,&
55 : transfer_pw2rs,&
56 : transfer_rs2pw
57 : USE shell_potential_types, ONLY: shell_kind_type
58 : USE structure_factor_types, ONLY: structure_factor_type
59 : USE structure_factors, ONLY: structure_factor_allocate,&
60 : structure_factor_deallocate,&
61 : structure_factor_init
62 : #include "./base/base_uses.f90"
63 :
64 : IMPLICIT NONE
65 :
66 : PRIVATE
67 : PUBLIC :: pme_evaluate
68 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pme'
69 :
70 : CONTAINS
71 :
72 : ! **************************************************************************************************
73 : !> \brief ...
74 : !> \param ewald_env ...
75 : !> \param ewald_pw ...
76 : !> \param box ...
77 : !> \param particle_set ...
78 : !> \param vg_coulomb ...
79 : !> \param fg_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 (15-Mar-2001) : New electrostatic calculation and pressure tensor
90 : !> JGH (21-Mar-2001) : Complete rewrite
91 : !> JGH (21-Mar-2001) : Introduced real space density type for future
92 : !> parallelisation
93 : !> \author CJM and APSI
94 : ! **************************************************************************************************
95 1818 : SUBROUTINE pme_evaluate(ewald_env, ewald_pw, box, particle_set, vg_coulomb, &
96 1818 : fg_coulomb, pv_g, shell_particle_set, core_particle_set, &
97 3636 : fgshell_coulomb, fgcore_coulomb, use_virial, charges, atprop)
98 : TYPE(ewald_environment_type), POINTER :: ewald_env
99 : TYPE(ewald_pw_type), POINTER :: ewald_pw
100 : TYPE(cell_type), POINTER :: box
101 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
102 : REAL(KIND=dp), INTENT(OUT) :: vg_coulomb
103 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), &
104 : OPTIONAL :: fg_coulomb, pv_g
105 : TYPE(particle_type), DIMENSION(:), OPTIONAL, &
106 : POINTER :: shell_particle_set, core_particle_set
107 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), &
108 : OPTIONAL :: fgshell_coulomb, fgcore_coulomb
109 : LOGICAL, INTENT(IN) :: use_virial
110 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: charges
111 : TYPE(atprop_type), POINTER :: atprop
112 :
113 : CHARACTER(LEN=*), PARAMETER :: routineN = 'pme_evaluate'
114 :
115 : INTEGER :: handle, i, ig, ipart, j, nd(3), npart, &
116 : nshell, p1, p2
117 : LOGICAL :: is1_core, is2_core
118 : REAL(KIND=dp) :: alpha, dvols, fat1, ffa, ffb
119 : REAL(KIND=dp), DIMENSION(3) :: fat
120 : REAL(KIND=dp), DIMENSION(3, 3) :: f_stress, h_stress
121 : TYPE(dg_rho0_type), POINTER :: dg_rho0
122 : TYPE(dg_type), POINTER :: dg
123 : TYPE(mp_comm_type) :: group
124 : TYPE(pw_c1d_gs_type) :: phi_g, rhob_g
125 7272 : TYPE(pw_c1d_gs_type), DIMENSION(3) :: dphi_g
126 : TYPE(pw_grid_type), POINTER :: grid_b, grid_s
127 : TYPE(pw_poisson_type), POINTER :: poisson_env
128 : TYPE(pw_pool_type), POINTER :: pw_big_pool, pw_small_pool
129 : TYPE(pw_r3d_rs_type) :: phi_r, rhob_r, rhos1, rhos2
130 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
131 65448 : TYPE(realspace_grid_type), DIMENSION(3) :: drpot
132 : TYPE(realspace_grid_type), POINTER :: rden, rpot
133 : TYPE(structure_factor_type) :: exp_igr
134 :
135 1818 : CALL timeset(routineN, handle)
136 1818 : NULLIFY (poisson_env, rden)
137 1818 : CALL cite_reference(Darden1993)
138 1818 : CALL ewald_env_get(ewald_env, alpha=alpha, group=group)
139 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_big_pool, &
140 : pw_small_pool=pw_small_pool, rs_desc=rs_desc, &
141 1818 : poisson_env=poisson_env, dg=dg)
142 :
143 1818 : grid_b => pw_big_pool%pw_grid
144 1818 : grid_s => pw_small_pool%pw_grid
145 :
146 1818 : CALL dg_get(dg, dg_rho0=dg_rho0)
147 :
148 1818 : npart = SIZE(particle_set)
149 :
150 1818 : CALL structure_factor_init(exp_igr)
151 :
152 1818 : IF (PRESENT(shell_particle_set)) THEN
153 0 : CPASSERT(ASSOCIATED(shell_particle_set))
154 0 : CPASSERT(ASSOCIATED(core_particle_set))
155 0 : nshell = SIZE(shell_particle_set)
156 : CALL structure_factor_allocate(grid_s%bounds, npart, exp_igr, &
157 : allocate_centre=.TRUE., allocate_shell_e=.TRUE., &
158 0 : allocate_shell_centre=.TRUE., nshell=nshell)
159 :
160 : ELSE
161 : CALL structure_factor_allocate(grid_s%bounds, npart, exp_igr, &
162 1818 : allocate_centre=.TRUE.)
163 : END IF
164 :
165 1818 : CALL pw_small_pool%create_pw(rhos1)
166 1818 : CALL pw_small_pool%create_pw(rhos2)
167 :
168 38178 : ALLOCATE (rden)
169 1818 : CALL rs_grid_create(rden, rs_desc)
170 1818 : CALL rs_grid_set_box(grid_b, rs=rden)
171 1818 : CALL rs_grid_zero(rden)
172 :
173 1818 : CPASSERT(ASSOCIATED(box))
174 :
175 1818 : IF (rden%desc%parallel .AND. rden%desc%distributed) THEN
176 0 : CALL get_center(particle_set, box, exp_igr%centre, exp_igr%delta, grid_b%npts, 1)
177 : END IF
178 1818 : IF (PRESENT(shell_particle_set) .AND. rden%desc%parallel .AND. rden%desc%distributed) THEN
179 0 : CALL get_center(shell_particle_set, box, exp_igr%shell_centre, exp_igr%shell_delta, grid_b%npts, 1)
180 0 : CALL get_center(core_particle_set, box, exp_igr%core_centre, exp_igr%core_delta, grid_b%npts, 1)
181 : END IF
182 :
183 : !-------------- DENSITY CALCULATION ----------------
184 :
185 1818 : ipart = 0
186 29400 : DO
187 :
188 14700 : CALL set_list(particle_set, npart, exp_igr%centre, p1, rden, ipart, exp_igr%core_centre)
189 14700 : CALL set_list(particle_set, npart, exp_igr%centre, p2, rden, ipart, exp_igr%core_centre)
190 14700 : IF (p1 == 0 .AND. p2 == 0) EXIT
191 :
192 12882 : is1_core = (particle_set(p1)%shell_index /= 0)
193 12882 : IF (p2 /= 0) THEN
194 11757 : is2_core = (particle_set(p2)%shell_index /= 0)
195 : ELSE
196 1125 : is2_core = .FALSE.
197 : END IF
198 :
199 : ! calculate function on small boxes (we use double packing in FFT)
200 14700 : IF (is1_core .OR. is2_core) THEN
201 : CALL get_patch(dg, particle_set, exp_igr, box, p1, p2, grid_b, grid_s, &
202 : rhos1, rhos2, is1_core=is1_core, is2_core=is2_core, &
203 0 : core_particle_set=core_particle_set, charges=charges)
204 :
205 : ! add boxes to real space grid (big box)
206 0 : IF (is1_core) THEN
207 0 : CALL dg_sum_patch(rden, rhos1, exp_igr%core_centre(:, particle_set(p1)%shell_index))
208 : ELSE
209 0 : CALL dg_sum_patch(rden, rhos1, exp_igr%centre(:, p1))
210 : END IF
211 0 : IF (p2 /= 0 .AND. is2_core) THEN
212 0 : CALL dg_sum_patch(rden, rhos2, exp_igr%core_centre(:, particle_set(p2)%shell_index))
213 0 : ELSE IF (p2 /= 0) THEN
214 0 : CALL dg_sum_patch(rden, rhos2, exp_igr%centre(:, p2))
215 : END IF
216 : ELSE
217 : CALL get_patch(dg, particle_set, exp_igr, box, p1, p2, grid_b, grid_s, &
218 25562 : rhos1, rhos2, charges=charges)
219 : ! add boxes to real space grid (big box)
220 12882 : CALL dg_sum_patch(rden, rhos1, exp_igr%centre(:, p1))
221 12882 : IF (p2 /= 0) CALL dg_sum_patch(rden, rhos2, exp_igr%centre(:, p2))
222 : END IF
223 :
224 : END DO
225 1818 : IF (PRESENT(shell_particle_set)) THEN
226 0 : ipart = 0
227 0 : DO
228 0 : CALL set_list(shell_particle_set, nshell, exp_igr%shell_centre, p1, rpot, ipart)
229 0 : CALL set_list(shell_particle_set, nshell, exp_igr%shell_centre, p2, rpot, ipart)
230 0 : IF (p1 == 0 .AND. p2 == 0) EXIT
231 : ! calculate function on small boxes (we use double packing in FFT)
232 : CALL get_patch(dg, shell_particle_set, exp_igr, box, p1, p2, grid_b, grid_s, &
233 0 : rhos1, rhos2, is1_shell=.TRUE., is2_shell=.TRUE., charges=charges)
234 : ! add boxes to real space grid (big box)
235 0 : CALL dg_sum_patch(rpot, rhos1, exp_igr%shell_centre(:, p1))
236 0 : IF (p2 /= 0) CALL dg_sum_patch(rpot, rhos2, exp_igr%shell_centre(:, p2))
237 : END DO
238 : END IF
239 :
240 1818 : CALL pw_big_pool%create_pw(rhob_r)
241 1818 : CALL transfer_rs2pw(rden, rhob_r)
242 :
243 : !-------------- ELECTROSTATIC CALCULATION -----------
244 :
245 : ! allocate intermediate arrays
246 7272 : DO i = 1, 3
247 7272 : CALL pw_big_pool%create_pw(dphi_g(i))
248 : END DO
249 1818 : CALL pw_big_pool%create_pw(phi_r)
250 :
251 1818 : CALL pw_poisson_solve(poisson_env, rhob_r, vg_coulomb, phi_r, dphi_g, h_stress)
252 :
253 : ! atomic energies
254 1818 : IF (atprop%energy .OR. atprop%stress) THEN
255 202 : dvols = rhos1%pw_grid%dvol
256 4242 : ALLOCATE (rpot)
257 202 : CALL rs_grid_create(rpot, rs_desc)
258 202 : CALL transfer_pw2rs(rpot, phi_r)
259 202 : ipart = 0
260 808 : DO
261 404 : CALL set_list(particle_set, npart, exp_igr%centre, p1, rden, ipart, exp_igr%core_centre)
262 404 : CALL set_list(particle_set, npart, exp_igr%centre, p2, rden, ipart, exp_igr%core_centre)
263 404 : IF (p1 == 0 .AND. p2 == 0) EXIT
264 : ! integrate box and potential
265 : CALL get_patch(dg, particle_set, exp_igr, box, p1, p2, grid_b, grid_s, &
266 404 : rhos1, rhos2, charges=charges)
267 : ! add boxes to real space grid (big box)
268 202 : CALL dg_sum_patch_force_1d(rpot, rhos1, exp_igr%centre(:, p1), fat1)
269 202 : IF (atprop%energy) THEN
270 202 : atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols
271 : END IF
272 202 : IF (atprop%stress) THEN
273 202 : atprop%atstress(1, 1, p1) = atprop%atstress(1, 1, p1) + 0.5_dp*fat1*dvols
274 202 : atprop%atstress(2, 2, p1) = atprop%atstress(2, 2, p1) + 0.5_dp*fat1*dvols
275 202 : atprop%atstress(3, 3, p1) = atprop%atstress(3, 3, p1) + 0.5_dp*fat1*dvols
276 : END IF
277 404 : IF (p2 /= 0) THEN
278 101 : CALL dg_sum_patch_force_1d(rpot, rhos2, exp_igr%centre(:, p2), fat1)
279 101 : IF (atprop%energy) THEN
280 101 : atprop%atener(p2) = atprop%atener(p2) + 0.5_dp*fat1*dvols
281 : END IF
282 101 : IF (atprop%stress) THEN
283 101 : atprop%atstress(1, 1, p2) = atprop%atstress(1, 1, p2) + 0.5_dp*fat1*dvols
284 101 : atprop%atstress(2, 2, p2) = atprop%atstress(2, 2, p2) + 0.5_dp*fat1*dvols
285 101 : atprop%atstress(3, 3, p2) = atprop%atstress(3, 3, p2) + 0.5_dp*fat1*dvols
286 : END IF
287 : END IF
288 : END DO
289 202 : IF (atprop%stress) THEN
290 202 : CALL pw_big_pool%create_pw(phi_g)
291 202 : CALL pw_big_pool%create_pw(rhob_g)
292 202 : ffa = (0.5_dp/dg_rho0%zet(1))**2
293 202 : ffb = 1.0_dp/fourpi
294 808 : DO i = 1, 3
295 4325022 : DO ig = grid_b%first_gne0, grid_b%ngpts_cut_local
296 4324416 : phi_g%array(ig) = ffb*dphi_g(i)%array(ig)*(ffa*grid_b%gsq(ig) + 1.0_dp)
297 4325022 : phi_g%array(ig) = phi_g%array(ig)*poisson_env%green_fft%influence_fn%array(ig)
298 : END DO
299 606 : IF (grid_b%have_g0) phi_g%array(1) = 0.0_dp
300 2020 : DO j = 1, i
301 1212 : CALL pw_copy(phi_g, rhob_g)
302 1212 : nd = 0
303 1212 : nd(j) = 1
304 1212 : CALL pw_derive(rhob_g, nd)
305 1212 : CALL pw_transfer(rhob_g, rhob_r)
306 1212 : CALL transfer_pw2rs(rpot, rhob_r)
307 :
308 1212 : ipart = 0
309 606 : DO
310 2424 : CALL set_list(particle_set, npart, exp_igr%centre, p1, rden, ipart, exp_igr%core_centre)
311 2424 : CALL set_list(particle_set, npart, exp_igr%centre, p2, rden, ipart, exp_igr%core_centre)
312 2424 : IF (p1 == 0 .AND. p2 == 0) EXIT
313 : ! integrate box and potential
314 : CALL get_patch(dg, particle_set, exp_igr, box, p1, p2, grid_b, grid_s, &
315 2424 : rhos1, rhos2, charges=charges)
316 : ! add boxes to real space grid (big box)
317 1212 : CALL dg_sum_patch_force_1d(rpot, rhos1, exp_igr%centre(:, p1), fat1)
318 1212 : atprop%atstress(i, j, p1) = atprop%atstress(i, j, p1) + fat1*dvols
319 1212 : IF (i /= j) atprop%atstress(j, i, p1) = atprop%atstress(j, i, p1) + fat1*dvols
320 2424 : IF (p2 /= 0) THEN
321 606 : CALL dg_sum_patch_force_1d(rpot, rhos2, exp_igr%centre(:, p2), fat1)
322 606 : atprop%atstress(i, j, p2) = atprop%atstress(i, j, p2) + fat1*dvols
323 606 : IF (i /= j) atprop%atstress(j, i, p2) = atprop%atstress(j, i, p2) + fat1*dvols
324 : END IF
325 : END DO
326 : END DO
327 : END DO
328 202 : CALL pw_big_pool%give_back_pw(phi_g)
329 202 : CALL pw_big_pool%give_back_pw(rhob_g)
330 : END IF
331 202 : CALL rs_grid_release(rpot)
332 202 : DEALLOCATE (rpot)
333 : END IF
334 :
335 1818 : CALL pw_big_pool%give_back_pw(phi_r)
336 :
337 : !---------- END OF ELECTROSTATIC CALCULATION --------
338 :
339 : !------------- STRESS TENSOR CALCULATION ------------
340 :
341 1818 : IF ((use_virial) .AND. (PRESENT(pv_g))) THEN
342 808 : DO i = 1, 3
343 2020 : DO j = i, 3
344 1212 : f_stress(i, j) = pw_integral_a2b(dphi_g(i), dphi_g(j))
345 1818 : f_stress(j, i) = f_stress(i, j)
346 : END DO
347 : END DO
348 202 : ffa = (1.0_dp/fourpi)*(0.5_dp/dg_rho0%zet(1))**2
349 2626 : f_stress = -ffa*f_stress
350 2626 : pv_g = h_stress + f_stress
351 : END IF
352 :
353 : !--------END OF STRESS TENSOR CALCULATION -----------
354 :
355 7272 : DO i = 1, 3
356 5454 : CALL rs_grid_create(drpot(i), rs_desc)
357 5454 : CALL rs_grid_set_box(grid_b, rs=drpot(i))
358 5454 : CALL pw_transfer(dphi_g(i), rhob_r)
359 5454 : CALL pw_big_pool%give_back_pw(dphi_g(i))
360 7272 : CALL transfer_pw2rs(drpot(i), rhob_r)
361 : END DO
362 :
363 1818 : CALL pw_big_pool%give_back_pw(rhob_r)
364 :
365 : !----------------- FORCE CALCULATION ----------------
366 :
367 : ! initialize the forces
368 1818 : IF (PRESENT(fg_coulomb)) THEN
369 199026 : fg_coulomb = 0.0_dp
370 1818 : dvols = rhos1%pw_grid%dvol
371 :
372 1818 : ipart = 0
373 29400 : DO
374 :
375 14700 : CALL set_list(particle_set, npart, exp_igr%centre, p1, rden, ipart, exp_igr%core_centre)
376 14700 : CALL set_list(particle_set, npart, exp_igr%centre, p2, rden, ipart, exp_igr%core_centre)
377 14700 : IF (p1 == 0 .AND. p2 == 0) EXIT
378 :
379 12882 : is1_core = (particle_set(p1)%shell_index /= 0)
380 12882 : IF (p2 /= 0) THEN
381 11757 : is2_core = (particle_set(p2)%shell_index /= 0)
382 : ELSE
383 1125 : is2_core = .FALSE.
384 : END IF
385 :
386 : ! calculate function on small boxes (we use double packing in FFT)
387 :
388 : CALL get_patch_again(dg, particle_set, exp_igr, p1, p2, rhos1, rhos2, &
389 25562 : is1_core=is1_core, is2_core=is2_core, charges=charges)
390 :
391 : ! sum boxes on real space grids (big box)
392 12882 : IF (is1_core) THEN
393 : CALL dg_sum_patch_force_3d(drpot, rhos1, &
394 0 : exp_igr%core_centre(:, particle_set(p1)%shell_index), fat)
395 : fgcore_coulomb(1, particle_set(p1)%shell_index) = &
396 0 : fgcore_coulomb(1, particle_set(p1)%shell_index) - fat(1)*dvols
397 : fgcore_coulomb(2, particle_set(p1)%shell_index) = &
398 0 : fgcore_coulomb(2, particle_set(p1)%shell_index) - fat(2)*dvols
399 : fgcore_coulomb(3, particle_set(p1)%shell_index) = &
400 0 : fgcore_coulomb(3, particle_set(p1)%shell_index) - fat(3)*dvols
401 : ELSE
402 12882 : CALL dg_sum_patch_force_3d(drpot, rhos1, exp_igr%centre(:, p1), fat)
403 12882 : fg_coulomb(1, p1) = fg_coulomb(1, p1) - fat(1)*dvols
404 12882 : fg_coulomb(2, p1) = fg_coulomb(2, p1) - fat(2)*dvols
405 12882 : fg_coulomb(3, p1) = fg_coulomb(3, p1) - fat(3)*dvols
406 : END IF
407 14700 : IF (p2 /= 0 .AND. is2_core) THEN
408 : CALL dg_sum_patch_force_3d(drpot, rhos1, &
409 0 : exp_igr%core_centre(:, particle_set(p2)%shell_index), fat)
410 : fgcore_coulomb(1, particle_set(p2)%shell_index) = &
411 0 : fgcore_coulomb(1, particle_set(p2)%shell_index) - fat(1)*dvols
412 : fgcore_coulomb(2, particle_set(p2)%shell_index) = &
413 0 : fgcore_coulomb(2, particle_set(p2)%shell_index) - fat(2)*dvols
414 : fgcore_coulomb(3, particle_set(p2)%shell_index) = &
415 0 : fgcore_coulomb(3, particle_set(p2)%shell_index) - fat(3)*dvols
416 12882 : ELSEIF (p2 /= 0) THEN
417 11757 : CALL dg_sum_patch_force_3d(drpot, rhos2, exp_igr%centre(:, p2), fat)
418 11757 : fg_coulomb(1, p2) = fg_coulomb(1, p2) - fat(1)*dvols
419 11757 : fg_coulomb(2, p2) = fg_coulomb(2, p2) - fat(2)*dvols
420 11757 : fg_coulomb(3, p2) = fg_coulomb(3, p2) - fat(3)*dvols
421 : END IF
422 :
423 : END DO
424 : END IF
425 1818 : IF (PRESENT(fgshell_coulomb)) THEN
426 0 : fgshell_coulomb = 0.0_dp
427 0 : dvols = rhos1%pw_grid%dvol
428 :
429 0 : ipart = 0
430 0 : DO
431 0 : CALL set_list(shell_particle_set, nshell, exp_igr%shell_centre, p1, rden, ipart)
432 0 : CALL set_list(shell_particle_set, nshell, exp_igr%shell_centre, p2, rden, ipart)
433 0 : IF (p1 == 0 .AND. p2 == 0) EXIT
434 :
435 : ! calculate function on small boxes (we use double packing in FFT)
436 : CALL get_patch_again(dg, shell_particle_set, exp_igr, p1, p2, rhos1, rhos2, &
437 0 : is1_shell=.TRUE., is2_shell=.TRUE., charges=charges)
438 :
439 : ! sum boxes on real space grids (big box)
440 0 : CALL dg_sum_patch_force_3d(drpot, rhos1, exp_igr%shell_centre(:, p1), fat)
441 0 : fgshell_coulomb(1, p1) = fgshell_coulomb(1, p1) - fat(1)*dvols
442 0 : fgshell_coulomb(2, p1) = fgshell_coulomb(2, p1) - fat(2)*dvols
443 0 : fgshell_coulomb(3, p1) = fgshell_coulomb(3, p1) - fat(3)*dvols
444 0 : IF (p2 /= 0) THEN
445 0 : CALL dg_sum_patch_force_3d(drpot, rhos2, exp_igr%shell_centre(:, p2), fat)
446 0 : fgshell_coulomb(1, p2) = fgshell_coulomb(1, p2) - fat(1)*dvols
447 0 : fgshell_coulomb(2, p2) = fgshell_coulomb(2, p2) - fat(2)*dvols
448 0 : fgshell_coulomb(3, p2) = fgshell_coulomb(3, p2) - fat(3)*dvols
449 : END IF
450 : END DO
451 :
452 : END IF
453 : !--------------END OF FORCE CALCULATION -------------
454 :
455 : !------------------CLEANING UP ----------------------
456 :
457 1818 : CALL rs_grid_release(rden)
458 1818 : DEALLOCATE (rden)
459 7272 : DO i = 1, 3
460 7272 : CALL rs_grid_release(drpot(i))
461 : END DO
462 :
463 1818 : CALL pw_small_pool%give_back_pw(rhos1)
464 1818 : CALL pw_small_pool%give_back_pw(rhos2)
465 1818 : CALL structure_factor_deallocate(exp_igr)
466 :
467 1818 : CALL timestop(handle)
468 :
469 29088 : END SUBROUTINE pme_evaluate
470 :
471 : ! **************************************************************************************************
472 : !> \brief Calculates local density in a small box
473 : !> \param dg ...
474 : !> \param particle_set ...
475 : !> \param exp_igr ...
476 : !> \param box ...
477 : !> \param p1 ...
478 : !> \param p2 ...
479 : !> \param grid_b ...
480 : !> \param grid_s ...
481 : !> \param rhos1 ...
482 : !> \param rhos2 ...
483 : !> \param is1_core ...
484 : !> \param is2_core ...
485 : !> \param is1_shell ...
486 : !> \param is2_shell ...
487 : !> \param core_particle_set ...
488 : !> \param charges ...
489 : !> \par History
490 : !> JGH (23-Mar-2001) : Switch to integer from particle list pointers
491 : !> \author JGH (21-Mar-2001)
492 : ! **************************************************************************************************
493 14296 : SUBROUTINE get_patch(dg, particle_set, exp_igr, box, p1, p2, &
494 : grid_b, grid_s, rhos1, rhos2, is1_core, is2_core, is1_shell, &
495 14296 : is2_shell, core_particle_set, charges)
496 :
497 : TYPE(dg_type), POINTER :: dg
498 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
499 : TYPE(structure_factor_type) :: exp_igr
500 : TYPE(cell_type), POINTER :: box
501 : INTEGER, INTENT(IN) :: p1, p2
502 : TYPE(pw_grid_type), INTENT(IN) :: grid_b, grid_s
503 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rhos1, rhos2
504 : LOGICAL, OPTIONAL :: is1_core, is2_core, is1_shell, is2_shell
505 : TYPE(particle_type), DIMENSION(:), OPTIONAL, &
506 : POINTER :: core_particle_set
507 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: charges
508 :
509 14296 : COMPLEX(KIND=dp), DIMENSION(:), POINTER :: ex1, ex2, ey1, ey2, ez1, ez2
510 14296 : INTEGER, DIMENSION(:), POINTER :: center1, center2
511 : LOGICAL :: my_is1_core, my_is1_shell, my_is2_core, &
512 : my_is2_shell, use_charge_array
513 : REAL(KIND=dp) :: q1, q2
514 : REAL(KIND=dp), DIMENSION(3) :: r1, r2
515 : TYPE(atomic_kind_type), POINTER :: atomic_kind
516 : TYPE(dg_rho0_type), POINTER :: dg_rho0
517 : TYPE(pw_r3d_rs_type), POINTER :: rho0
518 : TYPE(shell_kind_type), POINTER :: shell
519 :
520 14296 : NULLIFY (shell)
521 14296 : use_charge_array = PRESENT(charges)
522 14296 : my_is1_core = .FALSE.
523 14296 : my_is2_core = .FALSE.
524 0 : IF (PRESENT(is1_core)) my_is1_core = is1_core
525 14296 : IF (PRESENT(is2_core)) my_is2_core = is2_core
526 14296 : IF (my_is1_core .OR. my_is2_core) THEN
527 0 : CPASSERT(PRESENT(core_particle_set))
528 : END IF
529 14296 : my_is1_shell = .FALSE.
530 14296 : my_is2_shell = .FALSE.
531 14296 : IF (PRESENT(is1_shell)) my_is1_shell = is1_shell
532 14296 : IF (PRESENT(is2_shell)) my_is2_shell = is2_shell
533 14296 : IF (my_is1_core .AND. my_is1_shell) THEN
534 0 : CPABORT("Shell-model: cannot be core and shell simultaneously")
535 : END IF
536 :
537 14296 : CALL dg_get(dg, dg_rho0=dg_rho0)
538 14296 : rho0 => dg_rho0%density
539 :
540 14296 : IF (my_is1_core) THEN
541 0 : r1 = core_particle_set(particle_set(p1)%shell_index)%r
542 : ELSE
543 57184 : r1 = particle_set(p1)%r
544 : END IF
545 14296 : atomic_kind => particle_set(p1)%atomic_kind
546 14296 : IF (my_is1_core) THEN
547 0 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
548 0 : q1 = shell%charge_core
549 14296 : ELSE IF (my_is1_shell) THEN
550 0 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
551 0 : q1 = shell%charge_shell
552 : ELSE
553 14296 : CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q1)
554 : END IF
555 14296 : IF (use_charge_array) q1 = charges(p1)
556 :
557 14296 : IF (my_is1_shell) THEN
558 0 : center1 => exp_igr%shell_centre(:, p1)
559 0 : ex1 => exp_igr%shell_ex(:, p1)
560 0 : ey1 => exp_igr%shell_ey(:, p1)
561 0 : ez1 => exp_igr%shell_ez(:, p1)
562 14296 : ELSEIF (my_is1_core) THEN
563 0 : center1 => exp_igr%core_centre(:, particle_set(p1)%shell_index)
564 0 : ex1 => exp_igr%core_ex(:, particle_set(p1)%shell_index)
565 0 : ey1 => exp_igr%core_ey(:, particle_set(p1)%shell_index)
566 0 : ez1 => exp_igr%core_ez(:, particle_set(p1)%shell_index)
567 : ELSE
568 14296 : center1 => exp_igr%centre(:, p1)
569 14296 : ex1 => exp_igr%ex(:, p1)
570 14296 : ey1 => exp_igr%ey(:, p1)
571 14296 : ez1 => exp_igr%ez(:, p1)
572 : END IF
573 :
574 14296 : CPASSERT(ASSOCIATED(box))
575 :
576 : CALL dg_get_strucfac(box%hmat, r1, grid_s%npts, grid_b%npts, center1, &
577 14296 : exp_igr%lb, ex1, ey1, ez1)
578 :
579 14296 : IF (p2 /= 0) THEN
580 12464 : IF (my_is2_core) THEN
581 0 : r2 = core_particle_set(particle_set(p2)%shell_index)%r
582 : ELSE
583 49856 : r2 = particle_set(p2)%r
584 : END IF
585 12464 : atomic_kind => particle_set(p2)%atomic_kind
586 12464 : IF (my_is2_core) THEN
587 0 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
588 0 : q2 = shell%charge_core
589 12464 : ELSE IF (my_is2_shell) THEN
590 0 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
591 0 : q2 = shell%charge_shell
592 : ELSE
593 12464 : CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q2)
594 : END IF
595 12464 : IF (use_charge_array) q2 = charges(p2)
596 :
597 12464 : IF (my_is2_shell) THEN
598 0 : center2 => exp_igr%shell_centre(:, p2)
599 0 : ex2 => exp_igr%shell_ex(:, p2)
600 0 : ey2 => exp_igr%shell_ey(:, p2)
601 0 : ez2 => exp_igr%shell_ez(:, p2)
602 12464 : ELSEIF (my_is2_core) THEN
603 0 : center2 => exp_igr%core_centre(:, particle_set(p2)%shell_index)
604 0 : ex2 => exp_igr%core_ex(:, particle_set(p2)%shell_index)
605 0 : ey2 => exp_igr%core_ey(:, particle_set(p2)%shell_index)
606 0 : ez2 => exp_igr%core_ez(:, particle_set(p2)%shell_index)
607 : ELSE
608 12464 : center2 => exp_igr%centre(:, p2)
609 12464 : ex2 => exp_igr%ex(:, p2)
610 12464 : ey2 => exp_igr%ey(:, p2)
611 12464 : ez2 => exp_igr%ez(:, p2)
612 : END IF
613 : CALL dg_get_strucfac(box%hmat, r2, grid_s%npts, grid_b%npts, center2, &
614 12464 : exp_igr%lb, ex2, ey2, ez2)
615 : END IF
616 :
617 14296 : IF (p2 == 0) THEN
618 1832 : CALL dg_get_patch(rho0, rhos1, q1, ex1, ey1, ez1)
619 : ELSE
620 12464 : CALL dg_get_patch(rho0, rhos1, rhos2, q1, q2, ex1, ey1, ez1, ex2, ey2, ez2)
621 : END IF
622 :
623 14296 : END SUBROUTINE get_patch
624 :
625 : ! **************************************************************************************************
626 : !> \brief ...
627 : !> \param dg ...
628 : !> \param particle_set ...
629 : !> \param exp_igr ...
630 : !> \param p1 ...
631 : !> \param p2 ...
632 : !> \param rhos1 ...
633 : !> \param rhos2 ...
634 : !> \param is1_core ...
635 : !> \param is2_core ...
636 : !> \param is1_shell ...
637 : !> \param is2_shell ...
638 : !> \param charges ...
639 : ! **************************************************************************************************
640 12882 : SUBROUTINE get_patch_again(dg, particle_set, exp_igr, p1, p2, rhos1, rhos2, is1_core, &
641 12882 : is2_core, is1_shell, is2_shell, charges)
642 :
643 : TYPE(dg_type), POINTER :: dg
644 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
645 : TYPE(structure_factor_type) :: exp_igr
646 : INTEGER, INTENT(IN) :: p1, p2
647 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rhos1, rhos2
648 : LOGICAL, OPTIONAL :: is1_core, is2_core, is1_shell, is2_shell
649 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: charges
650 :
651 12882 : COMPLEX(KIND=dp), DIMENSION(:), POINTER :: ex1, ex2, ey1, ey2, ez1, ez2
652 : LOGICAL :: my_is1_core, my_is1_shell, my_is2_core, &
653 : my_is2_shell, use_charge_array
654 : REAL(KIND=dp) :: q1, q2
655 : TYPE(atomic_kind_type), POINTER :: atomic_kind
656 : TYPE(dg_rho0_type), POINTER :: dg_rho0
657 : TYPE(pw_r3d_rs_type), POINTER :: rho0
658 : TYPE(shell_kind_type), POINTER :: shell
659 :
660 12882 : NULLIFY (shell)
661 12882 : use_charge_array = PRESENT(charges)
662 12882 : my_is1_core = .FALSE.
663 12882 : my_is2_core = .FALSE.
664 12882 : IF (PRESENT(is1_core)) my_is1_core = is1_core
665 12882 : IF (PRESENT(is2_core)) my_is2_core = is2_core
666 12882 : my_is1_shell = .FALSE.
667 12882 : my_is2_shell = .FALSE.
668 12882 : IF (PRESENT(is1_shell)) my_is1_shell = is1_shell
669 12882 : IF (PRESENT(is2_shell)) my_is2_shell = is2_shell
670 :
671 12882 : CALL dg_get(dg, dg_rho0=dg_rho0)
672 12882 : rho0 => dg_rho0%density
673 :
674 12882 : atomic_kind => particle_set(p1)%atomic_kind
675 12882 : IF (my_is1_core) THEN
676 0 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
677 0 : q1 = shell%charge_core
678 12882 : ELSE IF (my_is1_shell) THEN
679 0 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
680 0 : q1 = shell%charge_shell
681 : ELSE
682 12882 : CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q1)
683 : END IF
684 12882 : IF (use_charge_array) q1 = charges(p1)
685 12882 : IF (my_is1_core) THEN
686 0 : ex1 => exp_igr%core_ex(:, particle_set(p1)%shell_index)
687 0 : ey1 => exp_igr%core_ey(:, particle_set(p1)%shell_index)
688 0 : ez1 => exp_igr%core_ez(:, particle_set(p1)%shell_index)
689 12882 : ELSEIF (my_is1_shell) THEN
690 0 : ex1 => exp_igr%shell_ex(:, p1)
691 0 : ey1 => exp_igr%shell_ey(:, p1)
692 0 : ez1 => exp_igr%shell_ez(:, p1)
693 : ELSE
694 12882 : ex1 => exp_igr%ex(:, p1)
695 12882 : ey1 => exp_igr%ey(:, p1)
696 12882 : ez1 => exp_igr%ez(:, p1)
697 : END IF
698 :
699 12882 : IF (p2 /= 0) THEN
700 11757 : atomic_kind => particle_set(p2)%atomic_kind
701 11757 : IF (my_is2_core) THEN
702 0 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
703 0 : q2 = shell%charge_core
704 11757 : ELSE IF (my_is2_shell) THEN
705 0 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
706 0 : q2 = shell%charge_shell
707 : ELSE
708 11757 : CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q2)
709 : END IF
710 11757 : IF (use_charge_array) q2 = charges(p2)
711 11757 : IF (my_is2_core) THEN
712 0 : ex2 => exp_igr%core_ex(:, particle_set(p2)%shell_index)
713 0 : ey2 => exp_igr%core_ey(:, particle_set(p2)%shell_index)
714 0 : ez2 => exp_igr%core_ez(:, particle_set(p2)%shell_index)
715 11757 : ELSEIF (my_is2_shell) THEN
716 0 : ex2 => exp_igr%shell_ex(:, p2)
717 0 : ey2 => exp_igr%shell_ey(:, p2)
718 0 : ez2 => exp_igr%shell_ez(:, p2)
719 : ELSE
720 11757 : ex2 => exp_igr%ex(:, p2)
721 11757 : ey2 => exp_igr%ey(:, p2)
722 11757 : ez2 => exp_igr%ez(:, p2)
723 : END IF
724 : END IF
725 :
726 12882 : IF (p2 == 0) THEN
727 1125 : CALL dg_get_patch(rho0, rhos1, q1, ex1, ey1, ez1)
728 : ELSE
729 : CALL dg_get_patch(rho0, rhos1, rhos2, q1, q2, &
730 11757 : ex1, ey1, ez1, ex2, ey2, ez2)
731 : END IF
732 :
733 12882 : END SUBROUTINE get_patch_again
734 :
735 : END MODULE pme
736 :
|