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