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