Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \par Literature
10 : !> M. Krack, A. Gambirasio, and M. Parrinello,
11 : !> "Ab-initio x-ray scattering of liquid water",
12 : !> J. Chem. Phys. 117, 9409 (2002)
13 : !> \author Matthias Krack
14 : !> \date 30.11.2005
15 : ! **************************************************************************************************
16 : MODULE xray_diffraction
17 : USE atomic_kind_types, ONLY: atomic_kind_type
18 : USE bibliography, ONLY: Krack2002,&
19 : cite_reference
20 : USE cp_control_types, ONLY: dft_control_type
21 : USE gapw_gspace_reconstruction, ONLY: calculate_rhotot_elec_gspace
22 : USE kinds, ONLY: dp,&
23 : int_8
24 : USE mathconstants, ONLY: twopi
25 : USE memory_utilities, ONLY: reallocate
26 : USE message_passing, ONLY: mp_para_env_type
27 : USE particle_types, ONLY: particle_type
28 : USE physcon, ONLY: angstrom
29 : USE pw_env_types, ONLY: pw_env_get,&
30 : pw_env_type
31 : USE pw_grids, ONLY: get_pw_grid_info
32 : USE pw_methods, ONLY: pw_zero
33 : USE pw_pool_types, ONLY: pw_pool_type
34 : USE pw_types, ONLY: pw_c1d_gs_type
35 : USE qs_environment_types, ONLY: get_qs_env,&
36 : qs_environment_type
37 : USE qs_rho_atom_types, ONLY: rho_atom_type
38 : USE qs_rho_types, ONLY: qs_rho_type
39 : USE util, ONLY: sort
40 : #include "./base/base_uses.f90"
41 :
42 : IMPLICIT NONE
43 :
44 : PRIVATE
45 :
46 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xray_diffraction'
47 :
48 : PUBLIC :: calculate_rhotot_elec_gspace, &
49 : xray_diffraction_spectrum
50 :
51 : CONTAINS
52 :
53 : ! **************************************************************************************************
54 : !> \brief Calculate the coherent X-ray diffraction spectrum using the total
55 : !> electronic density in reciprocal space (g-space).
56 : !> \param qs_env ...
57 : !> \param unit_number ...
58 : !> \param q_max ...
59 : !> \date 30.11.2005
60 : !> \author Matthias Krack
61 : ! **************************************************************************************************
62 120 : SUBROUTINE xray_diffraction_spectrum(qs_env, unit_number, q_max)
63 :
64 : TYPE(qs_environment_type), POINTER :: qs_env
65 : INTEGER, INTENT(IN) :: unit_number
66 : REAL(KIND=dp), INTENT(IN) :: q_max
67 :
68 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xray_diffraction_spectrum'
69 : INTEGER, PARAMETER :: nblock = 100
70 :
71 : INTEGER :: handle, i, ig, ig_shell, ipe, ishell, &
72 : jg, ng, npe, nshell, nshell_gather
73 : INTEGER(KIND=int_8) :: ngpts
74 : INTEGER, DIMENSION(3) :: npts
75 30 : INTEGER, DIMENSION(:), POINTER :: aux_index, ng_shell, ng_shell_gather, &
76 30 : nshell_pe, offset_pe
77 : REAL(KIND=dp) :: cutoff, f, f2, q, rho_hard, rho_soft, &
78 : rho_total
79 : REAL(KIND=dp), DIMENSION(3) :: dg, dr
80 30 : REAL(KIND=dp), DIMENSION(:), POINTER :: f2sum, f2sum_gather, f4sum, f4sum_gather, fmax, &
81 30 : fmax_gather, fmin, fmin_gather, fsum, fsum_gather, gsq, q_shell, q_shell_gather
82 30 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
83 : TYPE(dft_control_type), POINTER :: dft_control
84 : TYPE(mp_para_env_type), POINTER :: para_env
85 30 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
86 : TYPE(pw_c1d_gs_type) :: rhotot_elec_gspace
87 : TYPE(pw_env_type), POINTER :: pw_env
88 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
89 : TYPE(qs_rho_type), POINTER :: rho
90 30 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
91 :
92 0 : CPASSERT(ASSOCIATED(qs_env))
93 :
94 30 : CALL timeset(routineN, handle)
95 :
96 30 : NULLIFY (atomic_kind_set)
97 30 : NULLIFY (aux_index)
98 30 : NULLIFY (auxbas_pw_pool)
99 30 : NULLIFY (dft_control)
100 30 : NULLIFY (f2sum)
101 30 : NULLIFY (f2sum_gather)
102 30 : NULLIFY (f4sum)
103 30 : NULLIFY (f4sum_gather)
104 30 : NULLIFY (fmax)
105 30 : NULLIFY (fmax_gather)
106 30 : NULLIFY (fmin)
107 30 : NULLIFY (fmin_gather)
108 30 : NULLIFY (fsum)
109 30 : NULLIFY (fsum_gather)
110 30 : NULLIFY (gsq)
111 30 : NULLIFY (ng_shell)
112 30 : NULLIFY (ng_shell_gather)
113 30 : NULLIFY (nshell_pe)
114 30 : NULLIFY (offset_pe)
115 30 : NULLIFY (para_env)
116 30 : NULLIFY (particle_set)
117 30 : NULLIFY (pw_env)
118 30 : NULLIFY (q_shell)
119 30 : NULLIFY (q_shell_gather)
120 30 : NULLIFY (rho)
121 30 : NULLIFY (rho_atom_set)
122 :
123 30 : CALL cite_reference(Krack2002)
124 :
125 : CALL get_qs_env(qs_env=qs_env, &
126 : atomic_kind_set=atomic_kind_set, &
127 : dft_control=dft_control, &
128 : para_env=para_env, &
129 : particle_set=particle_set, &
130 : pw_env=pw_env, &
131 : rho=rho, &
132 30 : rho_atom_set=rho_atom_set)
133 :
134 : CALL pw_env_get(pw_env=pw_env, &
135 30 : auxbas_pw_pool=auxbas_pw_pool)
136 :
137 30 : npe = para_env%num_pe
138 :
139 : ! Plane waves grid to assemble the total electronic density
140 :
141 30 : CALL auxbas_pw_pool%create_pw(pw=rhotot_elec_gspace)
142 30 : CALL pw_zero(rhotot_elec_gspace)
143 :
144 : CALL get_pw_grid_info(pw_grid=rhotot_elec_gspace%pw_grid, &
145 : dr=dr, &
146 : npts=npts, &
147 : cutoff=cutoff, &
148 : ngpts=ngpts, &
149 30 : gsquare=gsq)
150 :
151 120 : dg(:) = twopi/(npts(:)*dr(:))
152 :
153 : ! Build the total electronic density in reciprocal space
154 :
155 : CALL calculate_rhotot_elec_gspace(qs_env=qs_env, &
156 : auxbas_pw_pool=auxbas_pw_pool, &
157 : rhotot_elec_gspace=rhotot_elec_gspace, &
158 : q_max=q_max, &
159 : rho_hard=rho_hard, &
160 30 : rho_soft=rho_soft)
161 :
162 30 : rho_total = rho_hard + rho_soft
163 :
164 : ! Calculate the coherent X-ray spectrum
165 :
166 : ! Now we have to gather the data from all processes, since each
167 : ! process has only worked his sub-grid
168 :
169 : ! Scan the g-vector shells
170 :
171 30 : CALL reallocate(q_shell, 1, nblock)
172 30 : CALL reallocate(ng_shell, 1, nblock)
173 :
174 30 : ng = SIZE(gsq)
175 :
176 30 : jg = 1
177 30 : nshell = 1
178 30 : q_shell(1) = SQRT(gsq(1))
179 30 : ng_shell(1) = 1
180 :
181 232603 : DO ig = 2, ng
182 232603 : CPASSERT(gsq(ig) >= gsq(jg))
183 232603 : IF (ABS(gsq(ig) - gsq(jg)) > 1.0E-12_dp) THEN
184 5359 : nshell = nshell + 1
185 5359 : IF (nshell > SIZE(q_shell)) THEN
186 44 : CALL reallocate(q_shell, 1, SIZE(q_shell) + nblock)
187 44 : CALL reallocate(ng_shell, 1, SIZE(ng_shell) + nblock)
188 : END IF
189 5359 : q = SQRT(gsq(ig))
190 5359 : IF (q > q_max) THEN
191 30 : nshell = nshell - 1
192 30 : EXIT
193 : END IF
194 5329 : q_shell(nshell) = q
195 5329 : ng_shell(nshell) = 1
196 5329 : jg = ig
197 : ELSE
198 227244 : ng_shell(nshell) = ng_shell(nshell) + 1
199 : END IF
200 : END DO
201 :
202 30 : CALL reallocate(q_shell, 1, nshell)
203 30 : CALL reallocate(ng_shell, 1, nshell)
204 30 : CALL reallocate(fmin, 1, nshell)
205 30 : CALL reallocate(fmax, 1, nshell)
206 30 : CALL reallocate(fsum, 1, nshell)
207 30 : CALL reallocate(f2sum, 1, nshell)
208 30 : CALL reallocate(f4sum, 1, nshell)
209 :
210 30 : ig = 0
211 5389 : DO ishell = 1, nshell
212 5359 : fmin(ishell) = HUGE(0.0_dp)
213 5359 : fmax(ishell) = 0.0_dp
214 5359 : fsum(ishell) = 0.0_dp
215 5359 : f2sum(ishell) = 0.0_dp
216 5359 : f4sum(ishell) = 0.0_dp
217 237962 : DO ig_shell = 1, ng_shell(ishell)
218 232603 : f = ABS(rhotot_elec_gspace%array(ig + ig_shell))
219 232603 : fmin(ishell) = MIN(fmin(ishell), f)
220 232603 : fmax(ishell) = MAX(fmax(ishell), f)
221 232603 : fsum(ishell) = fsum(ishell) + f
222 232603 : f2 = f*f
223 232603 : f2sum(ishell) = f2sum(ishell) + f2
224 237962 : f4sum(ishell) = f4sum(ishell) + f2*f2
225 : END DO
226 5389 : ig = ig + ng_shell(ishell)
227 : END DO
228 :
229 30 : CALL reallocate(nshell_pe, 0, npe - 1)
230 30 : CALL reallocate(offset_pe, 0, npe - 1)
231 :
232 : ! Root (source) process gathers the number of shell of each process
233 :
234 90 : CALL para_env%gather(nshell, nshell_pe)
235 :
236 : ! Only the root process which has to print the full spectrum has to
237 : ! allocate here the receive buffers with their real sizes
238 :
239 30 : IF (unit_number > 0) THEN
240 45 : nshell_gather = SUM(nshell_pe)
241 15 : offset_pe(0) = 0
242 30 : DO ipe = 1, npe - 1
243 30 : offset_pe(ipe) = offset_pe(ipe - 1) + nshell_pe(ipe - 1)
244 : END DO
245 : ELSE
246 15 : nshell_gather = 1 ! dummy value for the non-root processes
247 : END IF
248 :
249 30 : CALL reallocate(q_shell_gather, 1, nshell_gather)
250 30 : CALL reallocate(ng_shell_gather, 1, nshell_gather)
251 30 : CALL reallocate(fmin_gather, 1, nshell_gather)
252 30 : CALL reallocate(fmax_gather, 1, nshell_gather)
253 30 : CALL reallocate(fsum_gather, 1, nshell_gather)
254 30 : CALL reallocate(f2sum_gather, 1, nshell_gather)
255 30 : CALL reallocate(f4sum_gather, 1, nshell_gather)
256 :
257 10883 : CALL para_env%gatherv(q_shell, q_shell_gather, nshell_pe, offset_pe)
258 10883 : CALL para_env%gatherv(ng_shell, ng_shell_gather, nshell_pe, offset_pe)
259 10883 : CALL para_env%gatherv(fmax, fmax_gather, nshell_pe, offset_pe)
260 10883 : CALL para_env%gatherv(fmin, fmin_gather, nshell_pe, offset_pe)
261 10883 : CALL para_env%gatherv(fsum, fsum_gather, nshell_pe, offset_pe)
262 10883 : CALL para_env%gatherv(f2sum, f2sum_gather, nshell_pe, offset_pe)
263 10883 : CALL para_env%gatherv(f4sum, f4sum_gather, nshell_pe, offset_pe)
264 :
265 30 : IF (ASSOCIATED(offset_pe)) THEN
266 30 : DEALLOCATE (offset_pe)
267 : END IF
268 :
269 30 : IF (ASSOCIATED(nshell_pe)) THEN
270 30 : DEALLOCATE (nshell_pe)
271 : END IF
272 :
273 : ! Print X-ray diffraction spectrum (I/O node only)
274 :
275 30 : IF (unit_number > 0) THEN
276 :
277 15 : CALL reallocate(aux_index, 1, nshell_gather)
278 :
279 : ! Sort the gathered shells
280 :
281 15 : CALL sort(q_shell_gather, nshell_gather, aux_index)
282 :
283 : ! Allocate final arrays of sufficient size, i.e. nshell_gather
284 : ! is always greater or equal the final nshell value
285 :
286 15 : CALL reallocate(q_shell, 1, nshell_gather)
287 15 : CALL reallocate(ng_shell, 1, nshell_gather)
288 15 : CALL reallocate(fmin, 1, nshell_gather)
289 15 : CALL reallocate(fmax, 1, nshell_gather)
290 15 : CALL reallocate(fsum, 1, nshell_gather)
291 15 : CALL reallocate(f2sum, 1, nshell_gather)
292 15 : CALL reallocate(f4sum, 1, nshell_gather)
293 :
294 15 : jg = 1
295 15 : nshell = 1
296 15 : q_shell(1) = q_shell_gather(1)
297 15 : i = aux_index(1)
298 15 : ng_shell(1) = ng_shell_gather(i)
299 15 : fmin(1) = fmin_gather(i)
300 15 : fmax(1) = fmax_gather(i)
301 15 : fsum(1) = fsum_gather(i)
302 15 : f2sum(1) = f2sum_gather(i)
303 15 : f4sum(1) = f4sum_gather(i)
304 :
305 5359 : DO ig = 2, nshell_gather
306 5344 : i = aux_index(ig)
307 5359 : IF (ABS(q_shell_gather(ig) - q_shell_gather(jg)) > 1.0E-12_dp) THEN
308 2796 : nshell = nshell + 1
309 2796 : q_shell(nshell) = q_shell_gather(ig)
310 2796 : ng_shell(nshell) = ng_shell_gather(i)
311 2796 : fmin(nshell) = fmin_gather(i)
312 2796 : fmax(nshell) = fmax_gather(i)
313 2796 : fsum(nshell) = fsum_gather(i)
314 2796 : f2sum(nshell) = f2sum_gather(i)
315 2796 : f4sum(nshell) = f4sum_gather(i)
316 2796 : jg = ig
317 : ELSE
318 2548 : ng_shell(nshell) = ng_shell(nshell) + ng_shell_gather(i)
319 2548 : fmin(nshell) = MIN(fmin(nshell), fmin_gather(i))
320 2548 : fmax(nshell) = MAX(fmax(nshell), fmax_gather(i))
321 2548 : fsum(nshell) = fsum(nshell) + fsum_gather(i)
322 2548 : f2sum(nshell) = f2sum(nshell) + f2sum_gather(i)
323 2548 : f4sum(nshell) = f4sum(nshell) + f4sum_gather(i)
324 : END IF
325 : END DO
326 :
327 : ! The auxiliary index array is no longer needed now
328 :
329 15 : IF (ASSOCIATED(aux_index)) THEN
330 15 : DEALLOCATE (aux_index)
331 : END IF
332 :
333 : ! Allocate the final arrays for printing with their real size
334 :
335 15 : CALL reallocate(q_shell, 1, nshell)
336 15 : CALL reallocate(ng_shell, 1, nshell)
337 15 : CALL reallocate(fmin, 1, nshell)
338 15 : CALL reallocate(fmax, 1, nshell)
339 15 : CALL reallocate(fsum, 1, nshell)
340 15 : CALL reallocate(f2sum, 1, nshell)
341 15 : CALL reallocate(f4sum, 1, nshell)
342 :
343 : ! Write the X-ray diffraction spectrum to the specified file
344 :
345 : WRITE (UNIT=unit_number, FMT="(A)") &
346 15 : "#", &
347 15 : "# Coherent X-ray diffraction spectrum", &
348 30 : "#"
349 : WRITE (UNIT=unit_number, FMT="(A,1X,F20.10)") &
350 15 : "# Soft electronic charge (G-space) :", rho_soft, &
351 15 : "# Hard electronic charge (G-space) :", rho_hard, &
352 15 : "# Total electronic charge (G-space):", rho_total, &
353 15 : "# Density cutoff [Rydberg] :", 2.0_dp*cutoff, &
354 15 : "# q(min) [1/Angstrom] :", q_shell(2)/angstrom, &
355 15 : "# q(max) [1/Angstrom] :", q_shell(nshell)/angstrom, &
356 30 : "# q(max) [1/Angstrom] (requested) :", q_max/angstrom
357 : WRITE (UNIT=unit_number, FMT="(A,2X,I8)") &
358 15 : "# Number of g-vectors (grid points):", ngpts, &
359 30 : "# Number of g-vector shells :", nshell
360 : WRITE (UNIT=unit_number, FMT="(A,3(1X,I6))") &
361 15 : "# Grid size (a,b,c) :", npts(1:3)
362 : WRITE (UNIT=unit_number, FMT="(A,3F7.3)") &
363 60 : "# dg [1/Angstrom] :", dg(1:3)/angstrom, &
364 75 : "# dr [Angstrom] :", dr(1:3)*angstrom
365 : WRITE (UNIT=unit_number, FMT="(A)") &
366 15 : "#", &
367 : "# shell points q [1/A] <|F(q)|^2> Min(|F(q)|)"// &
368 30 : " Max(|F(q)|) <|F(q)|>^2 <|F(q)|^4>"
369 :
370 2826 : DO ishell = 1, nshell
371 : WRITE (UNIT=unit_number, FMT="(T2,I6,2X,I6,5(1X,F15.6),1X,ES15.6)") &
372 2811 : ishell, &
373 2811 : ng_shell(ishell), &
374 2811 : q_shell(ishell)/angstrom, &
375 2811 : f2sum(ishell)/REAL(ng_shell(ishell), KIND=dp), &
376 2811 : fmin(ishell), &
377 2811 : fmax(ishell), &
378 2811 : (fsum(ishell)/REAL(ng_shell(ishell), KIND=dp))**2, &
379 5637 : f4sum(ishell)/REAL(ng_shell(ishell), KIND=dp)
380 : END DO
381 :
382 : END IF
383 :
384 : ! Release work storage
385 :
386 30 : IF (ASSOCIATED(fmin)) THEN
387 30 : DEALLOCATE (fmin)
388 : END IF
389 :
390 30 : IF (ASSOCIATED(fmax)) THEN
391 30 : DEALLOCATE (fmax)
392 : END IF
393 :
394 30 : IF (ASSOCIATED(fsum)) THEN
395 30 : DEALLOCATE (fsum)
396 : END IF
397 :
398 30 : IF (ASSOCIATED(f2sum)) THEN
399 30 : DEALLOCATE (f2sum)
400 : END IF
401 :
402 30 : IF (ASSOCIATED(f4sum)) THEN
403 30 : DEALLOCATE (f4sum)
404 : END IF
405 :
406 30 : IF (ASSOCIATED(ng_shell)) THEN
407 30 : DEALLOCATE (ng_shell)
408 : END IF
409 :
410 30 : IF (ASSOCIATED(q_shell)) THEN
411 30 : DEALLOCATE (q_shell)
412 : END IF
413 :
414 30 : IF (ASSOCIATED(fmin_gather)) THEN
415 30 : DEALLOCATE (fmin_gather)
416 : END IF
417 :
418 30 : IF (ASSOCIATED(fmax_gather)) THEN
419 30 : DEALLOCATE (fmax_gather)
420 : END IF
421 :
422 30 : IF (ASSOCIATED(fsum_gather)) THEN
423 30 : DEALLOCATE (fsum_gather)
424 : END IF
425 :
426 30 : IF (ASSOCIATED(f2sum_gather)) THEN
427 30 : DEALLOCATE (f2sum_gather)
428 : END IF
429 :
430 30 : IF (ASSOCIATED(f4sum_gather)) THEN
431 30 : DEALLOCATE (f4sum_gather)
432 : END IF
433 :
434 30 : IF (ASSOCIATED(ng_shell_gather)) THEN
435 30 : DEALLOCATE (ng_shell_gather)
436 : END IF
437 :
438 30 : IF (ASSOCIATED(q_shell_gather)) THEN
439 30 : DEALLOCATE (q_shell_gather)
440 : END IF
441 :
442 30 : CALL auxbas_pw_pool%give_back_pw(rhotot_elec_gspace)
443 :
444 30 : CALL timestop(handle)
445 :
446 30 : END SUBROUTINE xray_diffraction_spectrum
447 :
448 : END MODULE xray_diffraction
|
