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 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 : get_atomic_kind
19 : USE basis_set_types, ONLY: get_gto_basis_set,&
20 : gto_basis_set_type
21 : USE bibliography, ONLY: Krack2002,&
22 : cite_reference
23 : USE cell_types, ONLY: cell_type,&
24 : pbc
25 : USE cp_control_types, ONLY: dft_control_type
26 : USE kinds, ONLY: dp,&
27 : int_8
28 : USE mathconstants, ONLY: pi,&
29 : twopi
30 : USE memory_utilities, ONLY: reallocate
31 : USE message_passing, ONLY: mp_para_env_type
32 : USE orbital_pointers, ONLY: indco,&
33 : nco,&
34 : ncoset,&
35 : nso,&
36 : nsoset
37 : USE orbital_transformation_matrices, ONLY: orbtramat
38 : USE particle_types, ONLY: particle_type
39 : USE paw_basis_types, ONLY: get_paw_basis_info
40 : USE physcon, ONLY: angstrom
41 : USE pw_env_types, ONLY: pw_env_get,&
42 : pw_env_type
43 : USE pw_grids, ONLY: get_pw_grid_info
44 : USE pw_methods, ONLY: pw_axpy,&
45 : pw_integrate_function,&
46 : pw_scale,&
47 : pw_transfer,&
48 : pw_zero
49 : USE pw_pool_types, ONLY: pw_pool_type
50 : USE pw_types, ONLY: pw_c1d_gs_type,&
51 : pw_r3d_rs_type
52 : USE qs_environment_types, ONLY: get_qs_env,&
53 : qs_environment_type
54 : USE qs_kind_types, ONLY: get_qs_kind,&
55 : qs_kind_type
56 : USE qs_rho_atom_types, ONLY: get_rho_atom,&
57 : rho_atom_coeff,&
58 : rho_atom_type
59 : USE qs_rho_types, ONLY: qs_rho_get,&
60 : qs_rho_type
61 : USE util, ONLY: sort
62 : #include "./base/base_uses.f90"
63 :
64 : IMPLICIT NONE
65 :
66 : PRIVATE
67 :
68 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xray_diffraction'
69 :
70 : PUBLIC :: calculate_rhotot_elec_gspace, &
71 : xray_diffraction_spectrum
72 :
73 : CONTAINS
74 :
75 : ! **************************************************************************************************
76 : !> \brief Calculate the coherent X-ray diffraction spectrum using the total
77 : !> electronic density in reciprocal space (g-space).
78 : !> \param qs_env ...
79 : !> \param unit_number ...
80 : !> \param q_max ...
81 : !> \date 30.11.2005
82 : !> \author Matthias Krack
83 : ! **************************************************************************************************
84 120 : SUBROUTINE xray_diffraction_spectrum(qs_env, unit_number, q_max)
85 :
86 : TYPE(qs_environment_type), POINTER :: qs_env
87 : INTEGER, INTENT(IN) :: unit_number
88 : REAL(KIND=dp), INTENT(IN) :: q_max
89 :
90 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xray_diffraction_spectrum'
91 : INTEGER, PARAMETER :: nblock = 100
92 :
93 : INTEGER :: handle, i, ig, ig_shell, ipe, ishell, &
94 : jg, ng, npe, nshell, nshell_gather
95 : INTEGER(KIND=int_8) :: ngpts
96 : INTEGER, DIMENSION(3) :: npts
97 30 : INTEGER, DIMENSION(:), POINTER :: aux_index, ng_shell, ng_shell_gather, &
98 30 : nshell_pe, offset_pe
99 : REAL(KIND=dp) :: cutoff, f, f2, q, rho_hard, rho_soft, &
100 : rho_total
101 : REAL(KIND=dp), DIMENSION(3) :: dg, dr
102 30 : REAL(KIND=dp), DIMENSION(:), POINTER :: f2sum, f2sum_gather, f4sum, f4sum_gather, fmax, &
103 30 : fmax_gather, fmin, fmin_gather, fsum, fsum_gather, gsq, q_shell, q_shell_gather
104 30 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
105 : TYPE(dft_control_type), POINTER :: dft_control
106 : TYPE(mp_para_env_type), POINTER :: para_env
107 30 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
108 : TYPE(pw_c1d_gs_type) :: rhotot_elec_gspace
109 : TYPE(pw_env_type), POINTER :: pw_env
110 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
111 : TYPE(qs_rho_type), POINTER :: rho
112 30 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
113 :
114 0 : CPASSERT(ASSOCIATED(qs_env))
115 :
116 30 : CALL timeset(routineN, handle)
117 :
118 30 : NULLIFY (atomic_kind_set)
119 30 : NULLIFY (aux_index)
120 30 : NULLIFY (auxbas_pw_pool)
121 30 : NULLIFY (dft_control)
122 30 : NULLIFY (f2sum)
123 30 : NULLIFY (f2sum_gather)
124 30 : NULLIFY (f4sum)
125 30 : NULLIFY (f4sum_gather)
126 30 : NULLIFY (fmax)
127 30 : NULLIFY (fmax_gather)
128 30 : NULLIFY (fmin)
129 30 : NULLIFY (fmin_gather)
130 30 : NULLIFY (fsum)
131 30 : NULLIFY (fsum_gather)
132 30 : NULLIFY (gsq)
133 30 : NULLIFY (ng_shell)
134 30 : NULLIFY (ng_shell_gather)
135 30 : NULLIFY (nshell_pe)
136 30 : NULLIFY (offset_pe)
137 30 : NULLIFY (para_env)
138 30 : NULLIFY (particle_set)
139 30 : NULLIFY (pw_env)
140 30 : NULLIFY (q_shell)
141 30 : NULLIFY (q_shell_gather)
142 30 : NULLIFY (rho)
143 30 : NULLIFY (rho_atom_set)
144 :
145 30 : CALL cite_reference(Krack2002)
146 :
147 : CALL get_qs_env(qs_env=qs_env, &
148 : atomic_kind_set=atomic_kind_set, &
149 : dft_control=dft_control, &
150 : para_env=para_env, &
151 : particle_set=particle_set, &
152 : pw_env=pw_env, &
153 : rho=rho, &
154 30 : rho_atom_set=rho_atom_set)
155 :
156 : CALL pw_env_get(pw_env=pw_env, &
157 30 : auxbas_pw_pool=auxbas_pw_pool)
158 :
159 30 : npe = para_env%num_pe
160 :
161 : ! Plane waves grid to assemble the total electronic density
162 :
163 30 : CALL auxbas_pw_pool%create_pw(pw=rhotot_elec_gspace)
164 30 : CALL pw_zero(rhotot_elec_gspace)
165 :
166 : CALL get_pw_grid_info(pw_grid=rhotot_elec_gspace%pw_grid, &
167 : dr=dr, &
168 : npts=npts, &
169 : cutoff=cutoff, &
170 : ngpts=ngpts, &
171 30 : gsquare=gsq)
172 :
173 120 : dg(:) = twopi/(npts(:)*dr(:))
174 :
175 : ! Build the total electronic density in reciprocal space
176 :
177 : CALL calculate_rhotot_elec_gspace(qs_env=qs_env, &
178 : auxbas_pw_pool=auxbas_pw_pool, &
179 : rhotot_elec_gspace=rhotot_elec_gspace, &
180 : q_max=q_max, &
181 : rho_hard=rho_hard, &
182 30 : rho_soft=rho_soft)
183 :
184 30 : rho_total = rho_hard + rho_soft
185 :
186 : ! Calculate the coherent X-ray spectrum
187 :
188 : ! Now we have to gather the data from all processes, since each
189 : ! process has only worked his sub-grid
190 :
191 : ! Scan the g-vector shells
192 :
193 30 : CALL reallocate(q_shell, 1, nblock)
194 30 : CALL reallocate(ng_shell, 1, nblock)
195 :
196 30 : ng = SIZE(gsq)
197 :
198 30 : jg = 1
199 30 : nshell = 1
200 30 : q_shell(1) = SQRT(gsq(1))
201 30 : ng_shell(1) = 1
202 :
203 232603 : DO ig = 2, ng
204 232603 : CPASSERT(gsq(ig) >= gsq(jg))
205 232603 : IF (ABS(gsq(ig) - gsq(jg)) > 1.0E-12_dp) THEN
206 5359 : nshell = nshell + 1
207 5359 : IF (nshell > SIZE(q_shell)) THEN
208 44 : CALL reallocate(q_shell, 1, SIZE(q_shell) + nblock)
209 44 : CALL reallocate(ng_shell, 1, SIZE(ng_shell) + nblock)
210 : END IF
211 5359 : q = SQRT(gsq(ig))
212 5359 : IF (q > q_max) THEN
213 30 : nshell = nshell - 1
214 30 : EXIT
215 : END IF
216 5329 : q_shell(nshell) = q
217 5329 : ng_shell(nshell) = 1
218 5329 : jg = ig
219 : ELSE
220 227244 : ng_shell(nshell) = ng_shell(nshell) + 1
221 : END IF
222 : END DO
223 :
224 30 : CALL reallocate(q_shell, 1, nshell)
225 30 : CALL reallocate(ng_shell, 1, nshell)
226 30 : CALL reallocate(fmin, 1, nshell)
227 30 : CALL reallocate(fmax, 1, nshell)
228 30 : CALL reallocate(fsum, 1, nshell)
229 30 : CALL reallocate(f2sum, 1, nshell)
230 30 : CALL reallocate(f4sum, 1, nshell)
231 :
232 30 : ig = 0
233 5389 : DO ishell = 1, nshell
234 5359 : fmin(ishell) = HUGE(0.0_dp)
235 5359 : fmax(ishell) = 0.0_dp
236 5359 : fsum(ishell) = 0.0_dp
237 5359 : f2sum(ishell) = 0.0_dp
238 5359 : f4sum(ishell) = 0.0_dp
239 237962 : DO ig_shell = 1, ng_shell(ishell)
240 232603 : f = ABS(rhotot_elec_gspace%array(ig + ig_shell))
241 232603 : fmin(ishell) = MIN(fmin(ishell), f)
242 232603 : fmax(ishell) = MAX(fmax(ishell), f)
243 232603 : fsum(ishell) = fsum(ishell) + f
244 232603 : f2 = f*f
245 232603 : f2sum(ishell) = f2sum(ishell) + f2
246 237962 : f4sum(ishell) = f4sum(ishell) + f2*f2
247 : END DO
248 5389 : ig = ig + ng_shell(ishell)
249 : END DO
250 :
251 30 : CALL reallocate(nshell_pe, 0, npe - 1)
252 30 : CALL reallocate(offset_pe, 0, npe - 1)
253 :
254 : ! Root (source) process gathers the number of shell of each process
255 :
256 90 : CALL para_env%gather(nshell, nshell_pe)
257 :
258 : ! Only the root process which has to print the full spectrum has to
259 : ! allocate here the receive buffers with their real sizes
260 :
261 30 : IF (unit_number > 0) THEN
262 45 : nshell_gather = SUM(nshell_pe)
263 15 : offset_pe(0) = 0
264 30 : DO ipe = 1, npe - 1
265 30 : offset_pe(ipe) = offset_pe(ipe - 1) + nshell_pe(ipe - 1)
266 : END DO
267 : ELSE
268 15 : nshell_gather = 1 ! dummy value for the non-root processes
269 : END IF
270 :
271 30 : CALL reallocate(q_shell_gather, 1, nshell_gather)
272 30 : CALL reallocate(ng_shell_gather, 1, nshell_gather)
273 30 : CALL reallocate(fmin_gather, 1, nshell_gather)
274 30 : CALL reallocate(fmax_gather, 1, nshell_gather)
275 30 : CALL reallocate(fsum_gather, 1, nshell_gather)
276 30 : CALL reallocate(f2sum_gather, 1, nshell_gather)
277 30 : CALL reallocate(f4sum_gather, 1, nshell_gather)
278 :
279 10883 : CALL para_env%gatherv(q_shell, q_shell_gather, nshell_pe, offset_pe)
280 10883 : CALL para_env%gatherv(ng_shell, ng_shell_gather, nshell_pe, offset_pe)
281 10883 : CALL para_env%gatherv(fmax, fmax_gather, nshell_pe, offset_pe)
282 10883 : CALL para_env%gatherv(fmin, fmin_gather, nshell_pe, offset_pe)
283 10883 : CALL para_env%gatherv(fsum, fsum_gather, nshell_pe, offset_pe)
284 10883 : CALL para_env%gatherv(f2sum, f2sum_gather, nshell_pe, offset_pe)
285 10883 : CALL para_env%gatherv(f4sum, f4sum_gather, nshell_pe, offset_pe)
286 :
287 30 : IF (ASSOCIATED(offset_pe)) THEN
288 30 : DEALLOCATE (offset_pe)
289 : END IF
290 :
291 30 : IF (ASSOCIATED(nshell_pe)) THEN
292 30 : DEALLOCATE (nshell_pe)
293 : END IF
294 :
295 : ! Print X-ray diffraction spectrum (I/O node only)
296 :
297 30 : IF (unit_number > 0) THEN
298 :
299 15 : CALL reallocate(aux_index, 1, nshell_gather)
300 :
301 : ! Sort the gathered shells
302 :
303 15 : CALL sort(q_shell_gather, nshell_gather, aux_index)
304 :
305 : ! Allocate final arrays of sufficient size, i.e. nshell_gather
306 : ! is always greater or equal the final nshell value
307 :
308 15 : CALL reallocate(q_shell, 1, nshell_gather)
309 15 : CALL reallocate(ng_shell, 1, nshell_gather)
310 15 : CALL reallocate(fmin, 1, nshell_gather)
311 15 : CALL reallocate(fmax, 1, nshell_gather)
312 15 : CALL reallocate(fsum, 1, nshell_gather)
313 15 : CALL reallocate(f2sum, 1, nshell_gather)
314 15 : CALL reallocate(f4sum, 1, nshell_gather)
315 :
316 15 : jg = 1
317 15 : nshell = 1
318 15 : q_shell(1) = q_shell_gather(1)
319 15 : i = aux_index(1)
320 15 : ng_shell(1) = ng_shell_gather(i)
321 15 : fmin(1) = fmin_gather(i)
322 15 : fmax(1) = fmax_gather(i)
323 15 : fsum(1) = fsum_gather(i)
324 15 : f2sum(1) = f2sum_gather(i)
325 15 : f4sum(1) = f4sum_gather(i)
326 :
327 5359 : DO ig = 2, nshell_gather
328 5344 : i = aux_index(ig)
329 5359 : IF (ABS(q_shell_gather(ig) - q_shell_gather(jg)) > 1.0E-12_dp) THEN
330 2796 : nshell = nshell + 1
331 2796 : q_shell(nshell) = q_shell_gather(ig)
332 2796 : ng_shell(nshell) = ng_shell_gather(i)
333 2796 : fmin(nshell) = fmin_gather(i)
334 2796 : fmax(nshell) = fmax_gather(i)
335 2796 : fsum(nshell) = fsum_gather(i)
336 2796 : f2sum(nshell) = f2sum_gather(i)
337 2796 : f4sum(nshell) = f4sum_gather(i)
338 2796 : jg = ig
339 : ELSE
340 2548 : ng_shell(nshell) = ng_shell(nshell) + ng_shell_gather(i)
341 2548 : fmin(nshell) = MIN(fmin(nshell), fmin_gather(i))
342 2548 : fmax(nshell) = MAX(fmax(nshell), fmax_gather(i))
343 2548 : fsum(nshell) = fsum(nshell) + fsum_gather(i)
344 2548 : f2sum(nshell) = f2sum(nshell) + f2sum_gather(i)
345 2548 : f4sum(nshell) = f4sum(nshell) + f4sum_gather(i)
346 : END IF
347 : END DO
348 :
349 : ! The auxiliary index array is no longer needed now
350 :
351 15 : IF (ASSOCIATED(aux_index)) THEN
352 15 : DEALLOCATE (aux_index)
353 : END IF
354 :
355 : ! Allocate the final arrays for printing with their real size
356 :
357 15 : CALL reallocate(q_shell, 1, nshell)
358 15 : CALL reallocate(ng_shell, 1, nshell)
359 15 : CALL reallocate(fmin, 1, nshell)
360 15 : CALL reallocate(fmax, 1, nshell)
361 15 : CALL reallocate(fsum, 1, nshell)
362 15 : CALL reallocate(f2sum, 1, nshell)
363 15 : CALL reallocate(f4sum, 1, nshell)
364 :
365 : ! Write the X-ray diffraction spectrum to the specified file
366 :
367 : WRITE (UNIT=unit_number, FMT="(A)") &
368 15 : "#", &
369 15 : "# Coherent X-ray diffraction spectrum", &
370 30 : "#"
371 : WRITE (UNIT=unit_number, FMT="(A,1X,F20.10)") &
372 15 : "# Soft electronic charge (G-space) :", rho_soft, &
373 15 : "# Hard electronic charge (G-space) :", rho_hard, &
374 15 : "# Total electronic charge (G-space):", rho_total, &
375 15 : "# Density cutoff [Rydberg] :", 2.0_dp*cutoff, &
376 15 : "# q(min) [1/Angstrom] :", q_shell(2)/angstrom, &
377 15 : "# q(max) [1/Angstrom] :", q_shell(nshell)/angstrom, &
378 30 : "# q(max) [1/Angstrom] (requested) :", q_max/angstrom
379 : WRITE (UNIT=unit_number, FMT="(A,2X,I8)") &
380 15 : "# Number of g-vectors (grid points):", ngpts, &
381 30 : "# Number of g-vector shells :", nshell
382 : WRITE (UNIT=unit_number, FMT="(A,3(1X,I6))") &
383 15 : "# Grid size (a,b,c) :", npts(1:3)
384 : WRITE (UNIT=unit_number, FMT="(A,3F7.3)") &
385 60 : "# dg [1/Angstrom] :", dg(1:3)/angstrom, &
386 75 : "# dr [Angstrom] :", dr(1:3)*angstrom
387 : WRITE (UNIT=unit_number, FMT="(A)") &
388 15 : "#", &
389 : "# shell points q [1/A] <|F(q)|^2> Min(|F(q)|)"// &
390 30 : " Max(|F(q)|) <|F(q)|>^2 <|F(q)|^4>"
391 :
392 2826 : DO ishell = 1, nshell
393 : WRITE (UNIT=unit_number, FMT="(T2,I6,2X,I6,5(1X,F15.6),1X,ES15.6)") &
394 2811 : ishell, &
395 2811 : ng_shell(ishell), &
396 2811 : q_shell(ishell)/angstrom, &
397 2811 : f2sum(ishell)/REAL(ng_shell(ishell), KIND=dp), &
398 2811 : fmin(ishell), &
399 2811 : fmax(ishell), &
400 2811 : (fsum(ishell)/REAL(ng_shell(ishell), KIND=dp))**2, &
401 5637 : f4sum(ishell)/REAL(ng_shell(ishell), KIND=dp)
402 : END DO
403 :
404 : END IF
405 :
406 : ! Release work storage
407 :
408 30 : IF (ASSOCIATED(fmin)) THEN
409 30 : DEALLOCATE (fmin)
410 : END IF
411 :
412 30 : IF (ASSOCIATED(fmax)) THEN
413 30 : DEALLOCATE (fmax)
414 : END IF
415 :
416 30 : IF (ASSOCIATED(fsum)) THEN
417 30 : DEALLOCATE (fsum)
418 : END IF
419 :
420 30 : IF (ASSOCIATED(f2sum)) THEN
421 30 : DEALLOCATE (f2sum)
422 : END IF
423 :
424 30 : IF (ASSOCIATED(f4sum)) THEN
425 30 : DEALLOCATE (f4sum)
426 : END IF
427 :
428 30 : IF (ASSOCIATED(ng_shell)) THEN
429 30 : DEALLOCATE (ng_shell)
430 : END IF
431 :
432 30 : IF (ASSOCIATED(q_shell)) THEN
433 30 : DEALLOCATE (q_shell)
434 : END IF
435 :
436 30 : IF (ASSOCIATED(fmin_gather)) THEN
437 30 : DEALLOCATE (fmin_gather)
438 : END IF
439 :
440 30 : IF (ASSOCIATED(fmax_gather)) THEN
441 30 : DEALLOCATE (fmax_gather)
442 : END IF
443 :
444 30 : IF (ASSOCIATED(fsum_gather)) THEN
445 30 : DEALLOCATE (fsum_gather)
446 : END IF
447 :
448 30 : IF (ASSOCIATED(f2sum_gather)) THEN
449 30 : DEALLOCATE (f2sum_gather)
450 : END IF
451 :
452 30 : IF (ASSOCIATED(f4sum_gather)) THEN
453 30 : DEALLOCATE (f4sum_gather)
454 : END IF
455 :
456 30 : IF (ASSOCIATED(ng_shell_gather)) THEN
457 30 : DEALLOCATE (ng_shell_gather)
458 : END IF
459 :
460 30 : IF (ASSOCIATED(q_shell_gather)) THEN
461 30 : DEALLOCATE (q_shell_gather)
462 : END IF
463 :
464 30 : CALL auxbas_pw_pool%give_back_pw(rhotot_elec_gspace)
465 :
466 30 : CALL timestop(handle)
467 :
468 30 : END SUBROUTINE xray_diffraction_spectrum
469 :
470 : ! **************************************************************************************************
471 : !> \brief The total electronic density in reciprocal space (g-space) is
472 : !> calculated.
473 : !> \param qs_env ...
474 : !> \param auxbas_pw_pool ...
475 : !> \param rhotot_elec_gspace ...
476 : !> \param q_max ...
477 : !> \param rho_hard ...
478 : !> \param rho_soft ...
479 : !> \param fsign ...
480 : !> \date 14.03.2008 (splitted from the routine xray_diffraction_spectrum)
481 : !> \author Matthias Krack
482 : !> \note This code assumes that the g-vectors are ordered (in gsq and %cc)
483 : ! **************************************************************************************************
484 108 : SUBROUTINE calculate_rhotot_elec_gspace(qs_env, auxbas_pw_pool, &
485 : rhotot_elec_gspace, q_max, rho_hard, &
486 : rho_soft, fsign)
487 :
488 : TYPE(qs_environment_type), POINTER :: qs_env
489 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
490 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rhotot_elec_gspace
491 : REAL(KIND=dp), INTENT(IN) :: q_max
492 : REAL(KIND=dp), INTENT(OUT) :: rho_hard, rho_soft
493 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: fsign
494 :
495 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_rhotot_elec_gspace'
496 :
497 : INTEGER :: atom, handle, iatom, ico, ico1_pgf, ico1_set, ikind, ipgf, iset, iso, iso1_pgf, &
498 : iso1_set, ison, ispin, jco, jco1_pgf, jco1_set, jpgf, jset, jso, jso1_pgf, jso1_set, &
499 : json, la, lb, maxco, maxso, na, natom, nb, ncoa, ncob, ncotot, nkind, nsatbas, nset, &
500 : nsoa, nsob, nsotot, nspin
501 36 : INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf, o2nindex
502 : LOGICAL :: orthorhombic, paw_atom
503 : REAL(KIND=dp) :: alpha, eps_rho_gspace, rho_total, scale, &
504 : volume
505 : REAL(KIND=dp), DIMENSION(3) :: ra
506 36 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: delta_cpc, pab, work, zet
507 36 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
508 : TYPE(cell_type), POINTER :: cell
509 : TYPE(dft_control_type), POINTER :: dft_control
510 : TYPE(gto_basis_set_type), POINTER :: basis_1c_set
511 36 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
512 : TYPE(pw_c1d_gs_type) :: rho_elec_gspace
513 36 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
514 36 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
515 : TYPE(qs_rho_type), POINTER :: rho
516 36 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: cpc_h, cpc_s
517 36 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
518 : TYPE(rho_atom_type), POINTER :: rho_atom
519 :
520 0 : CPASSERT(ASSOCIATED(qs_env))
521 36 : CPASSERT(ASSOCIATED(auxbas_pw_pool))
522 :
523 36 : CALL timeset(routineN, handle)
524 :
525 36 : NULLIFY (atom_list)
526 36 : NULLIFY (atomic_kind_set)
527 36 : NULLIFY (qs_kind_set)
528 36 : NULLIFY (cell)
529 36 : NULLIFY (cpc_h)
530 36 : NULLIFY (cpc_s)
531 36 : NULLIFY (delta_cpc)
532 36 : NULLIFY (dft_control)
533 36 : NULLIFY (lmax)
534 36 : NULLIFY (lmin)
535 36 : NULLIFY (npgf)
536 36 : NULLIFY (basis_1c_set)
537 36 : NULLIFY (pab)
538 36 : NULLIFY (particle_set)
539 36 : NULLIFY (rho, rho_r)
540 36 : NULLIFY (rho_atom)
541 36 : NULLIFY (rho_atom_set)
542 36 : NULLIFY (work)
543 36 : NULLIFY (zet)
544 :
545 : CALL get_qs_env(qs_env=qs_env, &
546 : atomic_kind_set=atomic_kind_set, &
547 : qs_kind_set=qs_kind_set, &
548 : cell=cell, &
549 : dft_control=dft_control, &
550 : particle_set=particle_set, &
551 : rho=rho, &
552 36 : rho_atom_set=rho_atom_set)
553 :
554 36 : CALL qs_rho_get(rho, rho_r=rho_r)
555 36 : eps_rho_gspace = dft_control%qs_control%eps_rho_gspace
556 36 : nkind = SIZE(atomic_kind_set)
557 36 : nspin = dft_control%nspins
558 :
559 : ! Load the soft contribution of the electronic density
560 :
561 36 : CALL auxbas_pw_pool%create_pw(pw=rho_elec_gspace)
562 :
563 36 : CALL pw_zero(rhotot_elec_gspace)
564 :
565 76 : DO ispin = 1, nspin
566 40 : CALL pw_zero(rho_elec_gspace)
567 40 : CALL pw_transfer(rho_r(ispin), rho_elec_gspace)
568 40 : IF (PRESENT(fsign) .AND. (ispin == 2)) THEN
569 2 : alpha = fsign
570 : ELSE
571 38 : alpha = 1.0_dp
572 : END IF
573 76 : CALL pw_axpy(rho_elec_gspace, rhotot_elec_gspace, alpha=alpha)
574 : END DO
575 :
576 : ! Release the auxiliary PW grid for the calculation of the soft
577 : ! contribution
578 :
579 36 : CALL auxbas_pw_pool%give_back_pw(rho_elec_gspace)
580 :
581 36 : rho_soft = pw_integrate_function(rhotot_elec_gspace, isign=-1)
582 :
583 : CALL get_pw_grid_info(pw_grid=rhotot_elec_gspace%pw_grid, &
584 : vol=volume, &
585 36 : orthorhombic=orthorhombic)
586 36 : IF (.NOT. orthorhombic) THEN
587 : CALL cp_abort(__LOCATION__, &
588 0 : "The calculation of XRD spectra for non-orthorhombic cells is not implemented")
589 : END IF
590 :
591 36 : CALL pw_scale(rhotot_elec_gspace, volume)
592 :
593 : ! Add the hard contribution of the electronic density
594 :
595 : ! Each process has to loop over all PAW atoms, since the g-space grid
596 : ! is already distributed over all processes
597 :
598 104 : DO ikind = 1, nkind
599 :
600 : CALL get_atomic_kind(atomic_kind_set(ikind), &
601 : atom_list=atom_list, &
602 68 : natom=natom)
603 :
604 : CALL get_qs_kind(qs_kind_set(ikind), &
605 : basis_set=basis_1c_set, &
606 : basis_type="GAPW_1C", &
607 68 : paw_atom=paw_atom)
608 :
609 68 : IF (.NOT. paw_atom) CYCLE ! no PAW atom: nothing to do
610 :
611 16 : CALL get_paw_basis_info(basis_1c_set, o2nindex=o2nindex, nsatbas=nsatbas)
612 :
613 : CALL get_gto_basis_set(gto_basis_set=basis_1c_set, &
614 : lmax=lmax, &
615 : lmin=lmin, &
616 : maxco=maxco, &
617 : maxso=maxso, &
618 : npgf=npgf, &
619 : nset=nset, &
620 16 : zet=zet)
621 :
622 16 : ncotot = maxco*nset
623 16 : nsotot = maxso*nset
624 16 : CALL reallocate(delta_cpc, 1, nsatbas, 1, nsatbas)
625 16 : CALL reallocate(pab, 1, ncotot, 1, ncotot)
626 16 : CALL reallocate(work, 1, maxso, 1, maxco)
627 :
628 40 : DO iatom = 1, natom
629 :
630 24 : atom = atom_list(iatom)
631 24 : rho_atom => rho_atom_set(atom)
632 :
633 : CALL get_rho_atom(rho_atom=rho_atom, &
634 : cpc_h=cpc_h, &
635 24 : cpc_s=cpc_s)
636 :
637 24 : ra(:) = pbc(particle_set(iatom)%r, cell)
638 :
639 12280 : delta_cpc = 0.0_dp
640 :
641 60 : DO ispin = 1, nspin
642 36 : IF (PRESENT(fsign) .AND. (ispin == 2)) THEN
643 6 : alpha = fsign
644 : ELSE
645 30 : alpha = 1.0_dp
646 : END IF
647 21388 : delta_cpc = delta_cpc + alpha*(cpc_h(ispin)%r_coef - cpc_s(ispin)%r_coef)
648 : END DO
649 :
650 24 : scale = 1.0_dp
651 :
652 120 : DO iset = 1, nset
653 80 : ico1_set = (iset - 1)*maxco + 1
654 80 : iso1_set = (iset - 1)*maxso + 1
655 80 : ncoa = ncoset(lmax(iset))
656 80 : nsoa = nsoset(lmax(iset))
657 392 : DO jset = 1, nset
658 288 : jco1_set = (jset - 1)*maxco + 1
659 288 : jso1_set = (jset - 1)*maxso + 1
660 288 : ncob = ncoset(lmax(jset))
661 288 : nsob = nsoset(lmax(jset))
662 1136 : DO ipgf = 1, npgf(iset)
663 768 : ico1_pgf = ico1_set + (ipgf - 1)*ncoa
664 768 : iso1_pgf = iso1_set + (ipgf - 1)*nsoa
665 3120 : DO jpgf = 1, npgf(jset)
666 2064 : jco1_pgf = jco1_set + (jpgf - 1)*ncob
667 2064 : jso1_pgf = jso1_set + (jpgf - 1)*nsob
668 2064 : ico = ico1_pgf + ncoset(lmin(iset) - 1)
669 2064 : iso = iso1_pgf + nsoset(lmin(iset) - 1)
670 :
671 : ! Transformation spherical to Cartesian
672 :
673 4832 : DO la = lmin(iset), lmax(iset)
674 2768 : jco = jco1_pgf + ncoset(lmin(jset) - 1)
675 2768 : jso = jso1_pgf + nsoset(lmin(jset) - 1)
676 6496 : DO lb = lmin(jset), lmax(jset)
677 3728 : ison = o2nindex(iso)
678 3728 : json = o2nindex(jso)
679 : CALL dgemm("N", "N", nso(la), nco(lb), nso(lb), 1.0_dp, &
680 : delta_cpc(ison:ison + nso(la) - 1, json), SIZE(delta_cpc, 1), &
681 : orbtramat(lb)%slm, nso(lb), 0.0_dp, work, &
682 3728 : maxso)
683 : CALL dgemm("T", "N", nco(la), nco(lb), nso(la), 1.0_dp, &
684 : orbtramat(la)%slm, nso(la), work, maxso, &
685 3728 : 0.0_dp, pab(ico:ico + nco(la) - 1, jco), SIZE(pab, 1))
686 3728 : jco = jco + nco(lb)
687 6496 : jso = jso + nso(lb)
688 : END DO ! next lb
689 2768 : ico = ico + nco(la)
690 4832 : iso = iso + nso(la)
691 : END DO ! next la
692 :
693 : ! Collocate current product of primitive Cartesian functions
694 :
695 2064 : na = ico1_pgf - 1
696 2064 : nb = jco1_pgf - 1
697 :
698 : CALL collocate_pgf_product_gspace( &
699 : la_max=lmax(iset), &
700 : zeta=zet(ipgf, iset), &
701 : la_min=lmin(iset), &
702 : lb_max=lmax(jset), &
703 : zetb=zet(jpgf, jset), &
704 : lb_min=lmin(jset), &
705 : ra=ra, &
706 : rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
707 : rab2=0.0_dp, &
708 : scale=scale, &
709 : pab=pab, &
710 : na=na, &
711 : nb=nb, &
712 : eps_rho_gspace=eps_rho_gspace, &
713 : gsq_max=q_max*q_max, &
714 2832 : pw=rhotot_elec_gspace)
715 :
716 : END DO ! next primitive Gaussian function "jpgf"
717 : END DO ! next primitive Gaussian function "ipgf"
718 : END DO ! next shell set "jset"
719 : END DO ! next shell set "iset"
720 : END DO ! next atom "iatom" of atomic kind "ikind"
721 136 : DEALLOCATE (o2nindex)
722 : END DO ! next atomic kind "ikind"
723 :
724 36 : rho_total = pw_integrate_function(rhotot_elec_gspace, isign=-1)/volume
725 :
726 36 : rho_hard = rho_total - rho_soft
727 :
728 : ! Release work storage
729 :
730 36 : IF (ASSOCIATED(delta_cpc)) THEN
731 8 : DEALLOCATE (delta_cpc)
732 : END IF
733 :
734 36 : IF (ASSOCIATED(work)) THEN
735 8 : DEALLOCATE (work)
736 : END IF
737 :
738 36 : IF (ASSOCIATED(pab)) THEN
739 8 : DEALLOCATE (pab)
740 : END IF
741 :
742 36 : CALL timestop(handle)
743 :
744 36 : END SUBROUTINE calculate_rhotot_elec_gspace
745 :
746 : ! **************************************************************************************************
747 : !> \brief low level collocation of primitive gaussian functions in g-space
748 : !> \param la_max ...
749 : !> \param zeta ...
750 : !> \param la_min ...
751 : !> \param lb_max ...
752 : !> \param zetb ...
753 : !> \param lb_min ...
754 : !> \param ra ...
755 : !> \param rab ...
756 : !> \param rab2 ...
757 : !> \param scale ...
758 : !> \param pab ...
759 : !> \param na ...
760 : !> \param nb ...
761 : !> \param eps_rho_gspace ...
762 : !> \param gsq_max ...
763 : !> \param pw ...
764 : ! **************************************************************************************************
765 2064 : SUBROUTINE collocate_pgf_product_gspace(la_max, zeta, la_min, &
766 : lb_max, zetb, lb_min, &
767 : ra, rab, rab2, scale, pab, na, nb, &
768 : eps_rho_gspace, gsq_max, pw)
769 :
770 : ! NOTE: this routine is much slower than the real-space version of collocate_pgf_product
771 :
772 : INTEGER, INTENT(IN) :: la_max
773 : REAL(dp), INTENT(IN) :: zeta
774 : INTEGER, INTENT(IN) :: la_min, lb_max
775 : REAL(dp), INTENT(IN) :: zetb
776 : INTEGER, INTENT(IN) :: lb_min
777 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra, rab
778 : REAL(dp), INTENT(IN) :: rab2, scale
779 : REAL(dp), DIMENSION(:, :), POINTER :: pab
780 : INTEGER, INTENT(IN) :: na, nb
781 : REAL(dp), INTENT(IN) :: eps_rho_gspace, gsq_max
782 : TYPE(pw_c1d_gs_type), INTENT(IN) :: pw
783 :
784 : CHARACTER(LEN=*), PARAMETER :: routineN = 'collocate_pgf_product_gspace'
785 :
786 : COMPLEX(dp), DIMENSION(3) :: phasefactor
787 2064 : COMPLEX(dp), DIMENSION(:), POINTER :: rag, rbg
788 2064 : COMPLEX(dp), DIMENSION(:, :, :, :), POINTER :: cubeaxis
789 : INTEGER :: ax, ay, az, bx, by, bz, handle, i, ico, &
790 : ig, ig2, jco, jg, kg, la, lb, &
791 : lb_cube_min, lb_grid, ub_cube_max, &
792 : ub_grid
793 : INTEGER, DIMENSION(3) :: lb_cube, ub_cube
794 : REAL(dp) :: f, fa, fb, pij, prefactor, rzetp, &
795 : twozetp, zetp
796 : REAL(dp), DIMENSION(3) :: dg, expfactor, fap, fbp, rap, rbp, rp
797 2064 : REAL(dp), DIMENSION(:), POINTER :: g
798 :
799 2064 : CALL timeset(routineN, handle)
800 :
801 8256 : dg(:) = twopi/(pw%pw_grid%npts(:)*pw%pw_grid%dr(:))
802 :
803 2064 : zetp = zeta + zetb
804 2064 : rzetp = 1.0_dp/zetp
805 2064 : f = zetb*rzetp
806 8256 : rap(:) = f*rab(:)
807 8256 : rbp(:) = rap(:) - rab(:)
808 8256 : rp(:) = ra(:) + rap(:)
809 2064 : twozetp = 2.0_dp*zetp
810 8256 : fap(:) = twozetp*rap(:)
811 8256 : fbp(:) = twozetp*rbp(:)
812 :
813 2064 : prefactor = scale*SQRT((pi*rzetp)**3)*EXP(-zeta*f*rab2)
814 8256 : phasefactor(:) = EXP(CMPLX(0.0_dp, -rp(:)*dg(:), KIND=dp))
815 8256 : expfactor(:) = EXP(-0.25*rzetp*dg(:)*dg(:))
816 :
817 8256 : lb_cube(:) = pw%pw_grid%bounds(1, :)
818 8256 : ub_cube(:) = pw%pw_grid%bounds(2, :)
819 :
820 8256 : lb_cube_min = MINVAL(lb_cube(:))
821 8256 : ub_cube_max = MAXVAL(ub_cube(:))
822 :
823 2064 : NULLIFY (cubeaxis, g, rag, rbg)
824 :
825 2064 : CALL reallocate(cubeaxis, lb_cube_min, ub_cube_max, 1, 3, 0, la_max, 0, lb_max)
826 2064 : CALL reallocate(g, lb_cube_min, ub_cube_max)
827 2064 : CALL reallocate(rag, lb_cube_min, ub_cube_max)
828 2064 : CALL reallocate(rbg, lb_cube_min, ub_cube_max)
829 :
830 2064 : lb_grid = LBOUND(pw%array, 1)
831 2064 : ub_grid = UBOUND(pw%array, 1)
832 :
833 8256 : DO i = 1, 3
834 :
835 340560 : DO ig = lb_cube(i), ub_cube(i)
836 334368 : ig2 = ig*ig
837 340560 : cubeaxis(ig, i, 0, 0) = expfactor(i)**ig2*phasefactor(i)**ig
838 : END DO
839 :
840 8256 : IF (la_max > 0) THEN
841 145200 : DO ig = lb_cube(i), ub_cube(i)
842 142560 : g(ig) = REAL(ig, dp)*dg(i)
843 142560 : rag(ig) = CMPLX(fap(i), -g(ig), KIND=dp)
844 145200 : cubeaxis(ig, i, 1, 0) = rag(ig)*cubeaxis(ig, i, 0, 0)
845 : END DO
846 3168 : DO la = 2, la_max
847 528 : fa = REAL(la - 1, dp)*twozetp
848 31680 : DO ig = lb_cube(i), ub_cube(i)
849 : cubeaxis(ig, i, la, 0) = rag(ig)*cubeaxis(ig, i, la - 1, 0) + &
850 29040 : fa*cubeaxis(ig, i, la - 2, 0)
851 : END DO
852 : END DO
853 2640 : IF (lb_max > 0) THEN
854 66000 : fa = twozetp
855 66000 : DO ig = lb_cube(i), ub_cube(i)
856 64800 : rbg(ig) = CMPLX(fbp(i), -g(ig), KIND=dp)
857 64800 : cubeaxis(ig, i, 0, 1) = rbg(ig)*cubeaxis(ig, i, 0, 0)
858 : cubeaxis(ig, i, 1, 1) = rbg(ig)*cubeaxis(ig, i, 1, 0) + &
859 66000 : fa*cubeaxis(ig, i, 0, 0)
860 : END DO
861 1440 : DO lb = 2, lb_max
862 240 : fb = REAL(lb - 1, dp)*twozetp
863 14400 : DO ig = lb_cube(i), ub_cube(i)
864 : cubeaxis(ig, i, 0, lb) = rbg(ig)*cubeaxis(ig, i, 0, lb - 1) + &
865 12960 : fb*cubeaxis(ig, i, 0, lb - 2)
866 : cubeaxis(ig, i, 1, lb) = rbg(ig)*cubeaxis(ig, i, 1, lb - 1) + &
867 : fb*cubeaxis(ig, i, 1, lb - 2) + &
868 13200 : fa*cubeaxis(ig, i, 0, lb - 1)
869 : END DO
870 : END DO
871 1440 : DO la = 2, la_max
872 240 : fa = REAL(la, dp)*twozetp
873 13200 : DO ig = lb_cube(i), ub_cube(i)
874 : cubeaxis(ig, i, la, 1) = rbg(ig)*cubeaxis(ig, i, la, 0) + &
875 13200 : fa*cubeaxis(ig, i, la - 1, 0)
876 : END DO
877 1488 : DO lb = 2, lb_max
878 48 : fb = REAL(lb - 1, dp)*twozetp
879 2880 : DO ig = lb_cube(i), ub_cube(i)
880 : cubeaxis(ig, i, la, lb) = rbg(ig)*cubeaxis(ig, i, la, lb - 1) + &
881 : fb*cubeaxis(ig, i, la, lb - 2) + &
882 2640 : fa*cubeaxis(ig, i, la - 1, lb - 1)
883 : END DO
884 : END DO
885 : END DO
886 : END IF
887 : ELSE
888 3552 : IF (lb_max > 0) THEN
889 79200 : DO ig = lb_cube(i), ub_cube(i)
890 77760 : g(ig) = REAL(ig, dp)*dg(i)
891 77760 : rbg(ig) = CMPLX(fbp(i), -g(ig), KIND=dp)
892 79200 : cubeaxis(ig, i, 0, 1) = rbg(ig)*cubeaxis(ig, i, 0, 0)
893 : END DO
894 1728 : DO lb = 2, lb_max
895 288 : fb = REAL(lb - 1, dp)*twozetp
896 17280 : DO ig = lb_cube(i), ub_cube(i)
897 : cubeaxis(ig, i, 0, lb) = rbg(ig)*cubeaxis(ig, i, 0, lb - 1) + &
898 15840 : fb*cubeaxis(ig, i, 0, lb - 2)
899 : END DO
900 : END DO
901 : END IF
902 : END IF
903 :
904 : END DO
905 :
906 5184 : DO la = 0, la_max
907 9936 : DO lb = 0, lb_max
908 4752 : IF (la + lb == 0) CYCLE
909 2688 : fa = (1.0_dp/twozetp)**(la + lb)
910 13872 : DO i = 1, 3
911 448272 : DO ig = lb_cube(i), ub_cube(i)
912 443520 : cubeaxis(ig, i, la, lb) = fa*cubeaxis(ig, i, la, lb)
913 : END DO
914 : END DO
915 : END DO
916 : END DO
917 :
918 : ! Add the current primitive Gaussian function product to grid
919 :
920 7120 : DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
921 :
922 5056 : ax = indco(1, ico)
923 5056 : ay = indco(2, ico)
924 5056 : az = indco(3, ico)
925 :
926 19792 : DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
927 :
928 12672 : pij = prefactor*pab(na + ico, nb + jco)
929 :
930 12672 : IF (ABS(pij) < eps_rho_gspace) CYCLE
931 :
932 4916 : bx = indco(1, jco)
933 4916 : by = indco(2, jco)
934 4916 : bz = indco(3, jco)
935 :
936 387056484 : DO i = lb_grid, ub_grid
937 387046512 : IF (pw%pw_grid%gsq(i) > gsq_max) CYCLE
938 346285458 : ig = pw%pw_grid%g_hat(1, i)
939 346285458 : jg = pw%pw_grid%g_hat(2, i)
940 346285458 : kg = pw%pw_grid%g_hat(3, i)
941 : pw%array(i) = pw%array(i) + pij*cubeaxis(ig, 1, ax, bx)* &
942 : cubeaxis(jg, 2, ay, by)* &
943 387059184 : cubeaxis(kg, 3, az, bz)
944 : END DO
945 :
946 : END DO
947 :
948 : END DO
949 :
950 2064 : DEALLOCATE (cubeaxis)
951 2064 : DEALLOCATE (g)
952 2064 : DEALLOCATE (rag)
953 2064 : DEALLOCATE (rbg)
954 :
955 2064 : CALL timestop(handle)
956 :
957 2064 : END SUBROUTINE collocate_pgf_product_gspace
958 :
959 : END MODULE xray_diffraction
|
