(git:ccc2433)
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 ! **************************************************************************************************
17  USE atomic_kind_types, ONLY: atomic_kind_type,&
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
38  USE particle_types, ONLY: particle_type
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, &
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  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 
959 END MODULE xray_diffraction
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
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...
Definition: bibliography.F:28
integer, save, public krack2002
Definition: bibliography.F:43
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.
Definition: mathconstants.F:16
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
Definition: pw_env_types.F:14
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
Definition: pw_env_types.F:113
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 ...
Definition: pw_pool_types.F:24
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.
Definition: qs_kind_types.F:23
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...
Definition: qs_rho_types.F:18
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...
Definition: qs_rho_types.F:229
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.