(git:374b731)
Loading...
Searching...
No Matches
xray_diffraction.F
Go to the documentation of this file.
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 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! **************************************************************************************************
21 USE bibliography, ONLY: krack2002,&
22 cite_reference
23 USE cell_types, ONLY: cell_type,&
24 pbc
26 USE kinds, ONLY: dp,&
27 int_8
28 USE mathconstants, ONLY: pi,&
29 twopi
32 USE orbital_pointers, ONLY: indco,&
33 nco,&
34 ncoset,&
35 nso,&
36 nsoset
40 USE physcon, ONLY: angstrom
41 USE pw_env_types, ONLY: pw_env_get,&
43 USE pw_grids, ONLY: get_pw_grid_info
44 USE pw_methods, ONLY: pw_axpy,&
46 pw_scale,&
50 USE pw_types, ONLY: pw_c1d_gs_type,&
54 USE qs_kind_types, ONLY: get_qs_kind,&
59 USE qs_rho_types, ONLY: qs_rho_get,&
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
72
73CONTAINS
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 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 INTEGER, DIMENSION(:), POINTER :: aux_index, ng_shell, ng_shell_gather, &
98 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 REAL(kind=dp), DIMENSION(:), POINTER :: f2sum, f2sum_gather, f4sum, f4sum_gather, fmax, &
103 fmax_gather, fmin, fmin_gather, fsum, fsum_gather, gsq, q_shell, q_shell_gather
104 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 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 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
113
114 cpassert(ASSOCIATED(qs_env))
115
116 CALL timeset(routinen, handle)
117
118 NULLIFY (atomic_kind_set)
119 NULLIFY (aux_index)
120 NULLIFY (auxbas_pw_pool)
121 NULLIFY (dft_control)
122 NULLIFY (f2sum)
123 NULLIFY (f2sum_gather)
124 NULLIFY (f4sum)
125 NULLIFY (f4sum_gather)
126 NULLIFY (fmax)
127 NULLIFY (fmax_gather)
128 NULLIFY (fmin)
129 NULLIFY (fmin_gather)
130 NULLIFY (fsum)
131 NULLIFY (fsum_gather)
132 NULLIFY (gsq)
133 NULLIFY (ng_shell)
134 NULLIFY (ng_shell_gather)
135 NULLIFY (nshell_pe)
136 NULLIFY (offset_pe)
137 NULLIFY (para_env)
138 NULLIFY (particle_set)
139 NULLIFY (pw_env)
140 NULLIFY (q_shell)
141 NULLIFY (q_shell_gather)
142 NULLIFY (rho)
143 NULLIFY (rho_atom_set)
144
145 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 rho_atom_set=rho_atom_set)
155
156 CALL pw_env_get(pw_env=pw_env, &
157 auxbas_pw_pool=auxbas_pw_pool)
158
159 npe = para_env%num_pe
160
161 ! Plane waves grid to assemble the total electronic density
162
163 CALL auxbas_pw_pool%create_pw(pw=rhotot_elec_gspace)
164 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 gsquare=gsq)
172
173 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 rho_soft=rho_soft)
183
184 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 CALL reallocate(q_shell, 1, nblock)
194 CALL reallocate(ng_shell, 1, nblock)
195
196 ng = SIZE(gsq)
197
198 jg = 1
199 nshell = 1
200 q_shell(1) = sqrt(gsq(1))
201 ng_shell(1) = 1
202
203 DO ig = 2, ng
204 cpassert(gsq(ig) >= gsq(jg))
205 IF (abs(gsq(ig) - gsq(jg)) > 1.0e-12_dp) THEN
206 nshell = nshell + 1
207 IF (nshell > SIZE(q_shell)) THEN
208 CALL reallocate(q_shell, 1, SIZE(q_shell) + nblock)
209 CALL reallocate(ng_shell, 1, SIZE(ng_shell) + nblock)
210 END IF
211 q = sqrt(gsq(ig))
212 IF (q > q_max) THEN
213 nshell = nshell - 1
214 EXIT
215 END IF
216 q_shell(nshell) = q
217 ng_shell(nshell) = 1
218 jg = ig
219 ELSE
220 ng_shell(nshell) = ng_shell(nshell) + 1
221 END IF
222 END DO
223
224 CALL reallocate(q_shell, 1, nshell)
225 CALL reallocate(ng_shell, 1, nshell)
226 CALL reallocate(fmin, 1, nshell)
227 CALL reallocate(fmax, 1, nshell)
228 CALL reallocate(fsum, 1, nshell)
229 CALL reallocate(f2sum, 1, nshell)
230 CALL reallocate(f4sum, 1, nshell)
231
232 ig = 0
233 DO ishell = 1, nshell
234 fmin(ishell) = huge(0.0_dp)
235 fmax(ishell) = 0.0_dp
236 fsum(ishell) = 0.0_dp
237 f2sum(ishell) = 0.0_dp
238 f4sum(ishell) = 0.0_dp
239 DO ig_shell = 1, ng_shell(ishell)
240 f = abs(rhotot_elec_gspace%array(ig + ig_shell))
241 fmin(ishell) = min(fmin(ishell), f)
242 fmax(ishell) = max(fmax(ishell), f)
243 fsum(ishell) = fsum(ishell) + f
244 f2 = f*f
245 f2sum(ishell) = f2sum(ishell) + f2
246 f4sum(ishell) = f4sum(ishell) + f2*f2
247 END DO
248 ig = ig + ng_shell(ishell)
249 END DO
250
251 CALL reallocate(nshell_pe, 0, npe - 1)
252 CALL reallocate(offset_pe, 0, npe - 1)
253
254 ! Root (source) process gathers the number of shell of each process
255
256 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 IF (unit_number > 0) THEN
262 nshell_gather = sum(nshell_pe)
263 offset_pe(0) = 0
264 DO ipe = 1, npe - 1
265 offset_pe(ipe) = offset_pe(ipe - 1) + nshell_pe(ipe - 1)
266 END DO
267 ELSE
268 nshell_gather = 1 ! dummy value for the non-root processes
269 END IF
270
271 CALL reallocate(q_shell_gather, 1, nshell_gather)
272 CALL reallocate(ng_shell_gather, 1, nshell_gather)
273 CALL reallocate(fmin_gather, 1, nshell_gather)
274 CALL reallocate(fmax_gather, 1, nshell_gather)
275 CALL reallocate(fsum_gather, 1, nshell_gather)
276 CALL reallocate(f2sum_gather, 1, nshell_gather)
277 CALL reallocate(f4sum_gather, 1, nshell_gather)
278
279 CALL para_env%gatherv(q_shell, q_shell_gather, nshell_pe, offset_pe)
280 CALL para_env%gatherv(ng_shell, ng_shell_gather, nshell_pe, offset_pe)
281 CALL para_env%gatherv(fmax, fmax_gather, nshell_pe, offset_pe)
282 CALL para_env%gatherv(fmin, fmin_gather, nshell_pe, offset_pe)
283 CALL para_env%gatherv(fsum, fsum_gather, nshell_pe, offset_pe)
284 CALL para_env%gatherv(f2sum, f2sum_gather, nshell_pe, offset_pe)
285 CALL para_env%gatherv(f4sum, f4sum_gather, nshell_pe, offset_pe)
286
287 IF (ASSOCIATED(offset_pe)) THEN
288 DEALLOCATE (offset_pe)
289 END IF
290
291 IF (ASSOCIATED(nshell_pe)) THEN
292 DEALLOCATE (nshell_pe)
293 END IF
294
295 ! Print X-ray diffraction spectrum (I/O node only)
296
297 IF (unit_number > 0) THEN
298
299 CALL reallocate(aux_index, 1, nshell_gather)
300
301 ! Sort the gathered shells
302
303 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 CALL reallocate(q_shell, 1, nshell_gather)
309 CALL reallocate(ng_shell, 1, nshell_gather)
310 CALL reallocate(fmin, 1, nshell_gather)
311 CALL reallocate(fmax, 1, nshell_gather)
312 CALL reallocate(fsum, 1, nshell_gather)
313 CALL reallocate(f2sum, 1, nshell_gather)
314 CALL reallocate(f4sum, 1, nshell_gather)
315
316 jg = 1
317 nshell = 1
318 q_shell(1) = q_shell_gather(1)
319 i = aux_index(1)
320 ng_shell(1) = ng_shell_gather(i)
321 fmin(1) = fmin_gather(i)
322 fmax(1) = fmax_gather(i)
323 fsum(1) = fsum_gather(i)
324 f2sum(1) = f2sum_gather(i)
325 f4sum(1) = f4sum_gather(i)
326
327 DO ig = 2, nshell_gather
328 i = aux_index(ig)
329 IF (abs(q_shell_gather(ig) - q_shell_gather(jg)) > 1.0e-12_dp) THEN
330 nshell = nshell + 1
331 q_shell(nshell) = q_shell_gather(ig)
332 ng_shell(nshell) = ng_shell_gather(i)
333 fmin(nshell) = fmin_gather(i)
334 fmax(nshell) = fmax_gather(i)
335 fsum(nshell) = fsum_gather(i)
336 f2sum(nshell) = f2sum_gather(i)
337 f4sum(nshell) = f4sum_gather(i)
338 jg = ig
339 ELSE
340 ng_shell(nshell) = ng_shell(nshell) + ng_shell_gather(i)
341 fmin(nshell) = min(fmin(nshell), fmin_gather(i))
342 fmax(nshell) = max(fmax(nshell), fmax_gather(i))
343 fsum(nshell) = fsum(nshell) + fsum_gather(i)
344 f2sum(nshell) = f2sum(nshell) + f2sum_gather(i)
345 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 IF (ASSOCIATED(aux_index)) THEN
352 DEALLOCATE (aux_index)
353 END IF
354
355 ! Allocate the final arrays for printing with their real size
356
357 CALL reallocate(q_shell, 1, nshell)
358 CALL reallocate(ng_shell, 1, nshell)
359 CALL reallocate(fmin, 1, nshell)
360 CALL reallocate(fmax, 1, nshell)
361 CALL reallocate(fsum, 1, nshell)
362 CALL reallocate(f2sum, 1, nshell)
363 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 "#", &
369 "# Coherent X-ray diffraction spectrum", &
370 "#"
371 WRITE (unit=unit_number, fmt="(A,1X,F20.10)") &
372 "# Soft electronic charge (G-space) :", rho_soft, &
373 "# Hard electronic charge (G-space) :", rho_hard, &
374 "# Total electronic charge (G-space):", rho_total, &
375 "# Density cutoff [Rydberg] :", 2.0_dp*cutoff, &
376 "# q(min) [1/Angstrom] :", q_shell(2)/angstrom, &
377 "# q(max) [1/Angstrom] :", q_shell(nshell)/angstrom, &
378 "# q(max) [1/Angstrom] (requested) :", q_max/angstrom
379 WRITE (unit=unit_number, fmt="(A,2X,I8)") &
380 "# Number of g-vectors (grid points):", ngpts, &
381 "# Number of g-vector shells :", nshell
382 WRITE (unit=unit_number, fmt="(A,3(1X,I6))") &
383 "# Grid size (a,b,c) :", npts(1:3)
384 WRITE (unit=unit_number, fmt="(A,3F7.3)") &
385 "# dg [1/Angstrom] :", dg(1:3)/angstrom, &
386 "# dr [Angstrom] :", dr(1:3)*angstrom
387 WRITE (unit=unit_number, fmt="(A)") &
388 "#", &
389 "# shell points q [1/A] <|F(q)|^2> Min(|F(q)|)"// &
390 " Max(|F(q)|) <|F(q)|>^2 <|F(q)|^4>"
391
392 DO ishell = 1, nshell
393 WRITE (unit=unit_number, fmt="(T2,I6,2X,I6,5(1X,F15.6),1X,ES15.6)") &
394 ishell, &
395 ng_shell(ishell), &
396 q_shell(ishell)/angstrom, &
397 f2sum(ishell)/real(ng_shell(ishell), kind=dp), &
398 fmin(ishell), &
399 fmax(ishell), &
400 (fsum(ishell)/real(ng_shell(ishell), kind=dp))**2, &
401 f4sum(ishell)/real(ng_shell(ishell), kind=dp)
402 END DO
403
404 END IF
405
406 ! Release work storage
407
408 IF (ASSOCIATED(fmin)) THEN
409 DEALLOCATE (fmin)
410 END IF
411
412 IF (ASSOCIATED(fmax)) THEN
413 DEALLOCATE (fmax)
414 END IF
415
416 IF (ASSOCIATED(fsum)) THEN
417 DEALLOCATE (fsum)
418 END IF
419
420 IF (ASSOCIATED(f2sum)) THEN
421 DEALLOCATE (f2sum)
422 END IF
423
424 IF (ASSOCIATED(f4sum)) THEN
425 DEALLOCATE (f4sum)
426 END IF
427
428 IF (ASSOCIATED(ng_shell)) THEN
429 DEALLOCATE (ng_shell)
430 END IF
431
432 IF (ASSOCIATED(q_shell)) THEN
433 DEALLOCATE (q_shell)
434 END IF
435
436 IF (ASSOCIATED(fmin_gather)) THEN
437 DEALLOCATE (fmin_gather)
438 END IF
439
440 IF (ASSOCIATED(fmax_gather)) THEN
441 DEALLOCATE (fmax_gather)
442 END IF
443
444 IF (ASSOCIATED(fsum_gather)) THEN
445 DEALLOCATE (fsum_gather)
446 END IF
447
448 IF (ASSOCIATED(f2sum_gather)) THEN
449 DEALLOCATE (f2sum_gather)
450 END IF
451
452 IF (ASSOCIATED(f4sum_gather)) THEN
453 DEALLOCATE (f4sum_gather)
454 END IF
455
456 IF (ASSOCIATED(ng_shell_gather)) THEN
457 DEALLOCATE (ng_shell_gather)
458 END IF
459
460 IF (ASSOCIATED(q_shell_gather)) THEN
461 DEALLOCATE (q_shell_gather)
462 END IF
463
464 CALL auxbas_pw_pool%give_back_pw(rhotot_elec_gspace)
465
466 CALL timestop(handle)
467
468 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 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: delta_cpc, pab, work, zet
507 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
512 TYPE(pw_c1d_gs_type) :: rho_elec_gspace
513 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
514 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
515 TYPE(qs_rho_type), POINTER :: rho
516 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: cpc_h, cpc_s
517 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
518 TYPE(rho_atom_type), POINTER :: rho_atom
519
520 cpassert(ASSOCIATED(qs_env))
521 cpassert(ASSOCIATED(auxbas_pw_pool))
522
523 CALL timeset(routinen, handle)
524
525 NULLIFY (atom_list)
526 NULLIFY (atomic_kind_set)
527 NULLIFY (qs_kind_set)
528 NULLIFY (cell)
529 NULLIFY (cpc_h)
530 NULLIFY (cpc_s)
531 NULLIFY (delta_cpc)
532 NULLIFY (dft_control)
533 NULLIFY (lmax)
534 NULLIFY (lmin)
535 NULLIFY (npgf)
536 NULLIFY (basis_1c_set)
537 NULLIFY (pab)
538 NULLIFY (particle_set)
539 NULLIFY (rho, rho_r)
540 NULLIFY (rho_atom)
541 NULLIFY (rho_atom_set)
542 NULLIFY (work)
543 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 rho_atom_set=rho_atom_set)
553
554 CALL qs_rho_get(rho, rho_r=rho_r)
555 eps_rho_gspace = dft_control%qs_control%eps_rho_gspace
556 nkind = SIZE(atomic_kind_set)
557 nspin = dft_control%nspins
558
559 ! Load the soft contribution of the electronic density
560
561 CALL auxbas_pw_pool%create_pw(pw=rho_elec_gspace)
562
563 CALL pw_zero(rhotot_elec_gspace)
564
565 DO ispin = 1, nspin
566 CALL pw_zero(rho_elec_gspace)
567 CALL pw_transfer(rho_r(ispin), rho_elec_gspace)
568 IF (PRESENT(fsign) .AND. (ispin == 2)) THEN
569 alpha = fsign
570 ELSE
571 alpha = 1.0_dp
572 END IF
573 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 CALL auxbas_pw_pool%give_back_pw(rho_elec_gspace)
580
581 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 orthorhombic=orthorhombic)
586 IF (.NOT. orthorhombic) THEN
587 CALL cp_abort(__location__, &
588 "The calculation of XRD spectra for non-orthorhombic cells is not implemented")
589 END IF
590
591 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 DO ikind = 1, nkind
599
600 CALL get_atomic_kind(atomic_kind_set(ikind), &
601 atom_list=atom_list, &
602 natom=natom)
603
604 CALL get_qs_kind(qs_kind_set(ikind), &
605 basis_set=basis_1c_set, &
606 basis_type="GAPW_1C", &
607 paw_atom=paw_atom)
608
609 IF (.NOT. paw_atom) cycle ! no PAW atom: nothing to do
610
611 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 zet=zet)
621
622 ncotot = maxco*nset
623 nsotot = maxso*nset
624 CALL reallocate(delta_cpc, 1, nsatbas, 1, nsatbas)
625 CALL reallocate(pab, 1, ncotot, 1, ncotot)
626 CALL reallocate(work, 1, maxso, 1, maxco)
627
628 DO iatom = 1, natom
629
630 atom = atom_list(iatom)
631 rho_atom => rho_atom_set(atom)
632
633 CALL get_rho_atom(rho_atom=rho_atom, &
634 cpc_h=cpc_h, &
635 cpc_s=cpc_s)
636
637 ra(:) = pbc(particle_set(iatom)%r, cell)
638
639 delta_cpc = 0.0_dp
640
641 DO ispin = 1, nspin
642 IF (PRESENT(fsign) .AND. (ispin == 2)) THEN
643 alpha = fsign
644 ELSE
645 alpha = 1.0_dp
646 END IF
647 delta_cpc = delta_cpc + alpha*(cpc_h(ispin)%r_coef - cpc_s(ispin)%r_coef)
648 END DO
649
650 scale = 1.0_dp
651
652 DO iset = 1, nset
653 ico1_set = (iset - 1)*maxco + 1
654 iso1_set = (iset - 1)*maxso + 1
655 ncoa = ncoset(lmax(iset))
656 nsoa = nsoset(lmax(iset))
657 DO jset = 1, nset
658 jco1_set = (jset - 1)*maxco + 1
659 jso1_set = (jset - 1)*maxso + 1
660 ncob = ncoset(lmax(jset))
661 nsob = nsoset(lmax(jset))
662 DO ipgf = 1, npgf(iset)
663 ico1_pgf = ico1_set + (ipgf - 1)*ncoa
664 iso1_pgf = iso1_set + (ipgf - 1)*nsoa
665 DO jpgf = 1, npgf(jset)
666 jco1_pgf = jco1_set + (jpgf - 1)*ncob
667 jso1_pgf = jso1_set + (jpgf - 1)*nsob
668 ico = ico1_pgf + ncoset(lmin(iset) - 1)
669 iso = iso1_pgf + nsoset(lmin(iset) - 1)
670
671 ! Transformation spherical to Cartesian
672
673 DO la = lmin(iset), lmax(iset)
674 jco = jco1_pgf + ncoset(lmin(jset) - 1)
675 jso = jso1_pgf + nsoset(lmin(jset) - 1)
676 DO lb = lmin(jset), lmax(jset)
677 ison = o2nindex(iso)
678 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 maxso)
683 CALL dgemm("T", "N", nco(la), nco(lb), nso(la), 1.0_dp, &
684 orbtramat(la)%slm, nso(la), work, maxso, &
685 0.0_dp, pab(ico:ico + nco(la) - 1, jco), SIZE(pab, 1))
686 jco = jco + nco(lb)
687 jso = jso + nso(lb)
688 END DO ! next lb
689 ico = ico + nco(la)
690 iso = iso + nso(la)
691 END DO ! next la
692
693 ! Collocate current product of primitive Cartesian functions
694
695 na = ico1_pgf - 1
696 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 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 DEALLOCATE (o2nindex)
722 END DO ! next atomic kind "ikind"
723
724 rho_total = pw_integrate_function(rhotot_elec_gspace, isign=-1)/volume
725
726 rho_hard = rho_total - rho_soft
727
728 ! Release work storage
729
730 IF (ASSOCIATED(delta_cpc)) THEN
731 DEALLOCATE (delta_cpc)
732 END IF
733
734 IF (ASSOCIATED(work)) THEN
735 DEALLOCATE (work)
736 END IF
737
738 IF (ASSOCIATED(pab)) THEN
739 DEALLOCATE (pab)
740 END IF
741
742 CALL timestop(handle)
743
744 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 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 COMPLEX(dp), DIMENSION(:), POINTER :: rag, rbg
788 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 REAL(dp), DIMENSION(:), POINTER :: g
798
799 CALL timeset(routinen, handle)
800
801 dg(:) = twopi/(pw%pw_grid%npts(:)*pw%pw_grid%dr(:))
802
803 zetp = zeta + zetb
804 rzetp = 1.0_dp/zetp
805 f = zetb*rzetp
806 rap(:) = f*rab(:)
807 rbp(:) = rap(:) - rab(:)
808 rp(:) = ra(:) + rap(:)
809 twozetp = 2.0_dp*zetp
810 fap(:) = twozetp*rap(:)
811 fbp(:) = twozetp*rbp(:)
812
813 prefactor = scale*sqrt((pi*rzetp)**3)*exp(-zeta*f*rab2)
814 phasefactor(:) = exp(cmplx(0.0_dp, -rp(:)*dg(:), kind=dp))
815 expfactor(:) = exp(-0.25*rzetp*dg(:)*dg(:))
816
817 lb_cube(:) = pw%pw_grid%bounds(1, :)
818 ub_cube(:) = pw%pw_grid%bounds(2, :)
819
820 lb_cube_min = minval(lb_cube(:))
821 ub_cube_max = maxval(ub_cube(:))
822
823 NULLIFY (cubeaxis, g, rag, rbg)
824
825 CALL reallocate(cubeaxis, lb_cube_min, ub_cube_max, 1, 3, 0, la_max, 0, lb_max)
826 CALL reallocate(g, lb_cube_min, ub_cube_max)
827 CALL reallocate(rag, lb_cube_min, ub_cube_max)
828 CALL reallocate(rbg, lb_cube_min, ub_cube_max)
829
830 lb_grid = lbound(pw%array, 1)
831 ub_grid = ubound(pw%array, 1)
832
833 DO i = 1, 3
834
835 DO ig = lb_cube(i), ub_cube(i)
836 ig2 = ig*ig
837 cubeaxis(ig, i, 0, 0) = expfactor(i)**ig2*phasefactor(i)**ig
838 END DO
839
840 IF (la_max > 0) THEN
841 DO ig = lb_cube(i), ub_cube(i)
842 g(ig) = real(ig, dp)*dg(i)
843 rag(ig) = cmplx(fap(i), -g(ig), kind=dp)
844 cubeaxis(ig, i, 1, 0) = rag(ig)*cubeaxis(ig, i, 0, 0)
845 END DO
846 DO la = 2, la_max
847 fa = real(la - 1, dp)*twozetp
848 DO ig = lb_cube(i), ub_cube(i)
849 cubeaxis(ig, i, la, 0) = rag(ig)*cubeaxis(ig, i, la - 1, 0) + &
850 fa*cubeaxis(ig, i, la - 2, 0)
851 END DO
852 END DO
853 IF (lb_max > 0) THEN
854 fa = twozetp
855 DO ig = lb_cube(i), ub_cube(i)
856 rbg(ig) = cmplx(fbp(i), -g(ig), kind=dp)
857 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 fa*cubeaxis(ig, i, 0, 0)
860 END DO
861 DO lb = 2, lb_max
862 fb = real(lb - 1, dp)*twozetp
863 DO ig = lb_cube(i), ub_cube(i)
864 cubeaxis(ig, i, 0, lb) = rbg(ig)*cubeaxis(ig, i, 0, lb - 1) + &
865 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 fa*cubeaxis(ig, i, 0, lb - 1)
869 END DO
870 END DO
871 DO la = 2, la_max
872 fa = real(la, dp)*twozetp
873 DO ig = lb_cube(i), ub_cube(i)
874 cubeaxis(ig, i, la, 1) = rbg(ig)*cubeaxis(ig, i, la, 0) + &
875 fa*cubeaxis(ig, i, la - 1, 0)
876 END DO
877 DO lb = 2, lb_max
878 fb = real(lb - 1, dp)*twozetp
879 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 fa*cubeaxis(ig, i, la - 1, lb - 1)
883 END DO
884 END DO
885 END DO
886 END IF
887 ELSE
888 IF (lb_max > 0) THEN
889 DO ig = lb_cube(i), ub_cube(i)
890 g(ig) = real(ig, dp)*dg(i)
891 rbg(ig) = cmplx(fbp(i), -g(ig), kind=dp)
892 cubeaxis(ig, i, 0, 1) = rbg(ig)*cubeaxis(ig, i, 0, 0)
893 END DO
894 DO lb = 2, lb_max
895 fb = real(lb - 1, dp)*twozetp
896 DO ig = lb_cube(i), ub_cube(i)
897 cubeaxis(ig, i, 0, lb) = rbg(ig)*cubeaxis(ig, i, 0, lb - 1) + &
898 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 DO la = 0, la_max
907 DO lb = 0, lb_max
908 IF (la + lb == 0) cycle
909 fa = (1.0_dp/twozetp)**(la + lb)
910 DO i = 1, 3
911 DO ig = lb_cube(i), ub_cube(i)
912 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 DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
921
922 ax = indco(1, ico)
923 ay = indco(2, ico)
924 az = indco(3, ico)
925
926 DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
927
928 pij = prefactor*pab(na + ico, nb + jco)
929
930 IF (abs(pij) < eps_rho_gspace) cycle
931
932 bx = indco(1, jco)
933 by = indco(2, jco)
934 bz = indco(3, jco)
935
936 DO i = lb_grid, ub_grid
937 IF (pw%pw_grid%gsq(i) > gsq_max) cycle
938 ig = pw%pw_grid%g_hat(1, i)
939 jg = pw%pw_grid%g_hat(2, i)
940 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 cubeaxis(kg, 3, az, bz)
944 END DO
945
946 END DO
947
948 END DO
949
950 DEALLOCATE (cubeaxis)
951 DEALLOCATE (g)
952 DEALLOCATE (rag)
953 DEALLOCATE (rbg)
954
955 CALL timestop(handle)
956
957 END SUBROUTINE collocate_pgf_product_gspace
958
959END MODULE xray_diffraction
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition atom.F:9
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius)
...
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public krack2002
Handles all functions related to the CELL.
Definition cell_types.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public int_8
Definition kinds.F:54
integer, parameter, public dp
Definition kinds.F:34
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
real(kind=dp), parameter, public twopi
Utility routines for the memory handling.
Interface to the message passing library MPI.
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public nco
integer, dimension(:), allocatable, public nsoset
integer, dimension(:), allocatable, public ncoset
integer, dimension(:, :), allocatable, public indco
integer, dimension(:), allocatable, public nso
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
type(orbtramat_type), dimension(:), pointer, public orbtramat
Define the data structure for the particle information.
subroutine, public get_paw_basis_info(basis_1c, o2nindex, n2oindex, nsatbas)
Return some info on the PAW basis derived from a GTO basis set.
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
This module defines the grid data type and some basic operations on it.
Definition pw_grids.F:36
subroutine, public get_pw_grid_info(pw_grid, id_nr, mode, vol, dvol, npts, ngpts, ngpts_cut, dr, cutoff, orthorhombic, gvectors, gsquare)
Access to information stored in the pw_grid_type.
Definition pw_grids.F:184
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_rho_atom(rho_atom, cpc_h, cpc_s, rho_rad_h, rho_rad_s, drho_rad_h, drho_rad_s, vrho_rad_h, vrho_rad_s, rho_rad_h_d, rho_rad_s_d, ga_vlocal_gb_h, ga_vlocal_gb_s, int_scr_h, int_scr_s)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
All kind of helpful little routines.
Definition util.F:14
subroutine, public xray_diffraction_spectrum(qs_env, unit_number, q_max)
Calculate the coherent X-ray diffraction spectrum using the total electronic density in reciprocal sp...
subroutine, public calculate_rhotot_elec_gspace(qs_env, auxbas_pw_pool, rhotot_elec_gspace, q_max, rho_hard, rho_soft, fsign)
The total electronic density in reciprocal space (g-space) is calculated.
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
stores all the informations relevant to an mpi environment
contained for different pw related things
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Provides all information about a quickstep kind.
keeps the density in various representations, keeping track of which ones are valid.