1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Utility methods to build 3-center integral tensors of various types. 8! ************************************************************************************************** 9MODULE qs_tensors 10 USE ai_contraction, ONLY: block_add 11 USE ai_contraction_sphi, ONLY: ab_contract,& 12 abc_contract 13 USE ai_overlap, ONLY: overlap_ab 14 USE ai_overlap3, ONLY: overlap3 15 USE atomic_kind_types, ONLY: atomic_kind_type 16 USE basis_set_types, ONLY: get_gto_basis_set,& 17 gto_basis_set_p_type 18 USE block_p_types, ONLY: block_p_type 19 USE cell_types, ONLY: cell_type 20 USE cp_control_types, ONLY: dft_control_type 21 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl 22 USE cp_files, ONLY: close_file,& 23 open_file 24 USE cp_para_types, ONLY: cp_para_env_type 25 USE dbcsr_api, ONLY: dbcsr_filter,& 26 dbcsr_finalize,& 27 dbcsr_get_block_p,& 28 dbcsr_has_symmetry,& 29 dbcsr_type,& 30 dbcsr_type_real_8 31 USE dbcsr_tensor_api, ONLY: & 32 dbcsr_t_blk_sizes, dbcsr_t_copy, dbcsr_t_create, dbcsr_t_destroy, & 33 dbcsr_t_distribution_destroy, dbcsr_t_distribution_new, dbcsr_t_distribution_type, & 34 dbcsr_t_filter, dbcsr_t_get_block, dbcsr_t_get_info, dbcsr_t_get_mapping_info, & 35 dbcsr_t_get_stored_coordinates, dbcsr_t_mp_environ_pgrid, dbcsr_t_pgrid_type, & 36 dbcsr_t_put_block, dbcsr_t_reserve_blocks, dbcsr_t_type, dbcsr_t_write_split_info 37 USE distribution_1d_types, ONLY: distribution_1d_type 38 USE distribution_2d_types, ONLY: distribution_2d_type 39 USE gamma, ONLY: init_md_ftable 40 USE input_constants, ONLY: do_potential_coulomb,& 41 do_potential_id,& 42 do_potential_short,& 43 do_potential_truncated 44 USE input_section_types, ONLY: section_vals_val_get 45 USE kinds, ONLY: default_string_length,& 46 dp 47 USE kpoint_types, ONLY: get_kpoint_info,& 48 kpoint_type 49 USE libint_2c_3c, ONLY: cutoff_screen_factor,& 50 eri_2center,& 51 eri_3center,& 52 libint_potential_type 53 USE libint_wrapper, ONLY: cp_libint_cleanup_2eri,& 54 cp_libint_cleanup_3eri,& 55 cp_libint_init_2eri,& 56 cp_libint_init_3eri,& 57 cp_libint_set_contrdepth,& 58 cp_libint_t 59 USE molecule_types, ONLY: molecule_type 60 USE orbital_pointers, ONLY: ncoset 61 USE particle_types, ONLY: particle_type 62 USE qs_environment_types, ONLY: get_qs_env,& 63 qs_environment_type 64 USE qs_kind_types, ONLY: qs_kind_type 65 USE qs_neighbor_list_types, ONLY: & 66 deallocate_neighbor_list_set, get_iterator_info, get_neighbor_list_set_p, & 67 neighbor_list_iterate, neighbor_list_iterator_create, neighbor_list_iterator_p_type, & 68 neighbor_list_iterator_release, neighbor_list_set_p_type, nl_sub_iterate 69 USE qs_neighbor_lists, ONLY: atom2d_build,& 70 atom2d_cleanup,& 71 build_neighbor_lists,& 72 local_atoms_type,& 73 pair_radius_setup 74 USE qs_tensors_types, ONLY: & 75 cyclic_tensor_dist, distribution_3d_destroy, distribution_3d_type, & 76 neighbor_list_3c_iterator_type, neighbor_list_3c_type, symmetric_ij, symmetric_ijk, & 77 symmetric_jk, symmetric_none, symmetrik_ik 78 USE t_c_g0, ONLY: get_lmax_init,& 79 init 80#include "./base/base_uses.f90" 81 82 IMPLICIT NONE 83 84 PRIVATE 85 86 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tensors' 87 88 PUBLIC :: build_3c_neighbor_lists, & 89 neighbor_list_3c_destroy, neighbor_list_2c_destroy, neighbor_list_3c_iterate, neighbor_list_3c_iterator_create, & 90 neighbor_list_3c_iterator_destroy, get_3c_iterator_info, build_3c_integrals, & 91 contiguous_tensor_dist, & 92 build_2c_neighbor_lists, build_2c_integrals, cutoff_screen_factor, tensor_change_pgrid 93 94 TYPE one_dim_int_array 95 INTEGER, DIMENSION(:), ALLOCATABLE :: array 96 END TYPE 97 98CONTAINS 99 100! ************************************************************************************************** 101!> \brief Build 2-center neighborlists adapted to different operators 102!> This mainly wraps build_neighbor_lists for consistency with build_3c_neighbor_lists 103!> \param ij_list 2c neighbor list for atom pairs i, j 104!> \param basis_i basis object for atoms i 105!> \param basis_j basis object for atoms j 106!> \param potential_parameter ... 107!> \param name name of 2c neighbor list 108!> \param qs_env ... 109!> \param sym_ij Symmetry in i, j (default .TRUE.) 110!> \param molecular ... 111!> \param dist_2d optionally a custom 2d distribution 112!> \param pot_to_rad which radius (1 for i, 2 for j) should be adapted to incorporate potential 113! ************************************************************************************************** 114 SUBROUTINE build_2c_neighbor_lists(ij_list, basis_i, basis_j, potential_parameter, name, qs_env, & 115 sym_ij, molecular, dist_2d, pot_to_rad) 116 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 117 POINTER :: ij_list 118 TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j 119 TYPE(libint_potential_type), INTENT(IN) :: potential_parameter 120 CHARACTER(LEN=*), INTENT(IN) :: name 121 TYPE(qs_environment_type), POINTER :: qs_env 122 LOGICAL, INTENT(IN), OPTIONAL :: sym_ij, molecular 123 TYPE(distribution_2d_type), OPTIONAL, POINTER :: dist_2d 124 INTEGER, INTENT(IN), OPTIONAL :: pot_to_rad 125 126 CHARACTER(len=*), PARAMETER :: routineN = 'build_2c_neighbor_lists', & 127 routineP = moduleN//':'//routineN 128 129 INTEGER :: ikind, nkind, pot_to_rad_prv 130 LOGICAL, ALLOCATABLE, DIMENSION(:) :: i_present, j_present 131 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius 132 REAL(kind=dp) :: subcells 133 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: i_radius, j_radius 134 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 135 TYPE(cell_type), POINTER :: cell 136 TYPE(distribution_1d_type), POINTER :: local_particles 137 TYPE(distribution_2d_type), POINTER :: dist_2d_prv 138 TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d 139 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 140 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 141 142 NULLIFY (atomic_kind_set, cell, local_particles, molecule_set, & 143 particle_set, dist_2d_prv) 144 145 IF (PRESENT(pot_to_rad)) THEN 146 pot_to_rad_prv = pot_to_rad 147 ELSE 148 pot_to_rad_prv = 1 149 ENDIF 150 151 CALL get_qs_env(qs_env, & 152 nkind=nkind, & 153 cell=cell, & 154 particle_set=particle_set, & 155 atomic_kind_set=atomic_kind_set, & 156 local_particles=local_particles, & 157 distribution_2d=dist_2d_prv, & 158 molecule_set=molecule_set) 159 160 CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells) 161 162 ALLOCATE (i_present(nkind), j_present(nkind)) 163 ALLOCATE (i_radius(nkind), j_radius(nkind)) 164 165 i_present = .FALSE. 166 j_present = .FALSE. 167 i_radius = 0.0_dp 168 j_radius = 0.0_dp 169 170 IF (PRESENT(dist_2d)) dist_2d_prv => dist_2d 171 172 ! Set up the radii, depending on the operator type 173 IF (potential_parameter%potential_type == do_potential_id) THEN 174 175 !overlap => use the kind radius for both i and j 176 DO ikind = 1, nkind 177 IF (ASSOCIATED(basis_i(ikind)%gto_basis_set)) THEN 178 i_present(ikind) = .TRUE. 179 CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=i_radius(ikind)) 180 END IF 181 IF (ASSOCIATED(basis_j(ikind)%gto_basis_set)) THEN 182 j_present(ikind) = .TRUE. 183 CALL get_gto_basis_set(basis_j(ikind)%gto_basis_set, kind_radius=j_radius(ikind)) 184 END IF 185 END DO 186 187 ELSE IF (potential_parameter%potential_type == do_potential_coulomb) THEN 188 189 !Coulomb operator, virtually infinite range => set j_radius to arbitrarily large number 190 DO ikind = 1, nkind 191 IF (ASSOCIATED(basis_i(ikind)%gto_basis_set)) THEN 192 i_present(ikind) = .TRUE. 193 IF (pot_to_rad_prv == 1) i_radius(ikind) = 1000000.0_dp 194 ENDIF 195 IF (ASSOCIATED(basis_j(ikind)%gto_basis_set)) THEN 196 j_present(ikind) = .TRUE. 197 IF (pot_to_rad_prv == 2) j_radius(ikind) = 1000000.0_dp 198 END IF 199 END DO !ikind 200 201 ELSE IF (potential_parameter%potential_type == do_potential_truncated .OR. & 202 potential_parameter%potential_type == do_potential_short) THEN 203 204 !Truncated coulomb/short range: set j_radius to r_cutoff + the kind_radii 205 DO ikind = 1, nkind 206 IF (ASSOCIATED(basis_i(ikind)%gto_basis_set)) THEN 207 i_present(ikind) = .TRUE. 208 CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=i_radius(ikind)) 209 IF (pot_to_rad_prv == 1) i_radius(ikind) = i_radius(ikind) + cutoff_screen_factor*potential_parameter%cutoff_radius 210 END IF 211 IF (ASSOCIATED(basis_j(ikind)%gto_basis_set)) THEN 212 j_present(ikind) = .TRUE. 213 CALL get_gto_basis_set(basis_j(ikind)%gto_basis_set, kind_radius=j_radius(ikind)) 214 IF (pot_to_rad_prv == 2) j_radius(ikind) = j_radius(ikind) + cutoff_screen_factor*potential_parameter%cutoff_radius 215 END IF 216 END DO 217 218 ELSE 219 CPABORT("Operator not implemented.") 220 END IF 221 222 ALLOCATE (pair_radius(nkind, nkind)) 223 pair_radius = 0.0_dp 224 CALL pair_radius_setup(i_present, j_present, i_radius, j_radius, pair_radius) 225 226 ALLOCATE (atom2d(nkind)) 227 228 CALL atom2d_build(atom2d, local_particles, dist_2d_prv, atomic_kind_set, & 229 molecule_set, molecule_only=.FALSE., particle_set=particle_set) 230 CALL build_neighbor_lists(ij_list, particle_set, atom2d, cell, pair_radius, subcells, & 231 symmetric=sym_ij, molecular=molecular, nlname=TRIM(name)) 232 233 CALL atom2d_cleanup(atom2d) 234 235 END SUBROUTINE 236 237! ************************************************************************************************** 238!> \brief ... 239!> \param ij_list ... 240! ************************************************************************************************** 241 SUBROUTINE neighbor_list_2c_destroy(ij_list) 242 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 243 POINTER :: ij_list 244 245 INTEGER :: i 246 247 DO i = 1, SIZE(ij_list) 248 CALL deallocate_neighbor_list_set(ij_list(i)%neighbor_list_set) 249 END DO 250 251 DEALLOCATE (ij_list) 252 253 END SUBROUTINE 254 255! ************************************************************************************************** 256!> \brief Build a 3-center neighbor list 257!> \param ijk_list 3c neighbor list for atom triples i, j, k 258!> \param basis_i basis object for atoms i 259!> \param basis_j basis object for atoms j 260!> \param basis_k basis object for atoms k 261!> \param dist_3d 3d distribution object 262!> \param potential_parameter ... 263!> \param name name of 3c neighbor list 264!> \param qs_env ... 265!> \param sym_ij Symmetry in i, j (default .FALSE.) 266!> \param sym_jk Symmetry in j, k (default .FALSE.) 267!> \param sym_ik Symmetry in i, k (default .FALSE.) 268!> \param molecular ??? not tested 269!> \param op_pos ... 270!> \param own_dist ... 271! ************************************************************************************************** 272 SUBROUTINE build_3c_neighbor_lists(ijk_list, basis_i, basis_j, basis_k, & 273 dist_3d, potential_parameter, name, qs_env, & 274 sym_ij, sym_jk, sym_ik, molecular, op_pos, own_dist) 275 TYPE(neighbor_list_3c_type), INTENT(OUT) :: ijk_list 276 TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j, basis_k 277 TYPE(distribution_3d_type), INTENT(IN) :: dist_3d 278 TYPE(libint_potential_type), INTENT(IN) :: potential_parameter 279 CHARACTER(LEN=*), INTENT(IN) :: name 280 TYPE(qs_environment_type), POINTER :: qs_env 281 LOGICAL, INTENT(IN), OPTIONAL :: sym_ij, sym_jk, sym_ik, molecular 282 INTEGER, INTENT(IN), OPTIONAL :: op_pos 283 LOGICAL, INTENT(IN), OPTIONAL :: own_dist 284 285 CHARACTER(len=*), PARAMETER :: routineN = 'build_3c_neighbor_lists', & 286 routineP = moduleN//':'//routineN 287 288 INTEGER :: handle, op_pos_prv, sym_level 289 TYPE(libint_potential_type) :: pot_par_1, pot_par_2 290 291 CALL timeset(routineN, handle) 292 293 IF (PRESENT(op_pos)) THEN 294 op_pos_prv = op_pos 295 ELSE 296 op_pos_prv = 1 297 ENDIF 298 299 SELECT CASE (op_pos_prv) 300 CASE (1) 301 pot_par_1 = potential_parameter 302 pot_par_2%potential_type = do_potential_id 303 CASE (2) 304 pot_par_2 = potential_parameter 305 pot_par_1%potential_type = do_potential_id 306 END SELECT 307 308 CALL build_2c_neighbor_lists(ijk_list%ij_list, basis_i, basis_j, pot_par_1, TRIM(name)//"_sub_1", & 309 qs_env, sym_ij=.FALSE., molecular=molecular, & 310 dist_2d=dist_3d%dist_2d_1, pot_to_rad=1) 311 312 CALL build_2c_neighbor_lists(ijk_list%jk_list, basis_j, basis_k, pot_par_2, TRIM(name)//"_sub_2", & 313 qs_env, sym_ij=.FALSE., molecular=molecular, & 314 dist_2d=dist_3d%dist_2d_2, pot_to_rad=2) 315 316 ijk_list%sym = symmetric_none 317 318 sym_level = 0 319 IF (PRESENT(sym_ij)) THEN 320 IF (sym_ij) THEN 321 ijk_list%sym = symmetric_ij 322 sym_level = sym_level + 1 323 ENDIF 324 ENDIF 325 326 IF (PRESENT(sym_jk)) THEN 327 IF (sym_jk) THEN 328 ijk_list%sym = symmetric_jk 329 sym_level = sym_level + 1 330 ENDIF 331 ENDIF 332 333 IF (PRESENT(sym_ik)) THEN 334 IF (sym_ik) THEN 335 ijk_list%sym = symmetrik_ik 336 sym_level = sym_level + 1 337 ENDIF 338 ENDIF 339 340 IF (sym_level >= 2) THEN 341 ijk_list%sym = symmetric_ijk 342 ENDIF 343 344 ijk_list%dist_3d = dist_3d 345 IF (PRESENT(own_dist)) THEN 346 ijk_list%owns_dist = own_dist 347 ELSE 348 ijk_list%owns_dist = .FALSE. 349 ENDIF 350 351 CALL timestop(handle) 352 END SUBROUTINE 353 354! ************************************************************************************************** 355!> \brief Symmetry criterion 356!> \param a ... 357!> \param b ... 358!> \return ... 359! ************************************************************************************************** 360 PURE FUNCTION include_symmetric(a, b) 361 INTEGER, INTENT(IN) :: a, b 362 LOGICAL :: include_symmetric 363 364 IF (a > b) THEN 365 include_symmetric = (MODULO(a + b, 2) /= 0) 366 ELSE 367 include_symmetric = (MODULO(a + b, 2) == 0) 368 END IF 369 370 END FUNCTION 371 372! ************************************************************************************************** 373!> \brief Destroy 3c neighborlist 374!> \param ijk_list ... 375! ************************************************************************************************** 376 SUBROUTINE neighbor_list_3c_destroy(ijk_list) 377 TYPE(neighbor_list_3c_type), INTENT(INOUT) :: ijk_list 378 379 CALL neighbor_list_2c_destroy(ijk_list%ij_list) 380 CALL neighbor_list_2c_destroy(ijk_list%jk_list) 381 382 IF (ijk_list%owns_dist) THEN 383 CALL distribution_3d_destroy(ijk_list%dist_3d) 384 ENDIF 385 386 END SUBROUTINE 387 388! ************************************************************************************************** 389!> \brief Create a 3-center neighborlist iterator 390!> \param iterator ... 391!> \param ijk_nl ... 392! ************************************************************************************************** 393 SUBROUTINE neighbor_list_3c_iterator_create(iterator, ijk_nl) 394 TYPE(neighbor_list_3c_iterator_type), INTENT(OUT) :: iterator 395 TYPE(neighbor_list_3c_type), INTENT(IN) :: ijk_nl 396 397 CHARACTER(len=*), PARAMETER :: routineN = 'neighbor_list_3c_iterator_create', & 398 routineP = moduleN//':'//routineN 399 400 INTEGER :: handle 401 402 CALL timeset(routineN, handle) 403 CALL neighbor_list_iterator_create(iterator%iter_ij, ijk_nl%ij_list) 404 CALL neighbor_list_iterator_create(iterator%iter_jk, ijk_nl%jk_list, search=.TRUE.) 405 iterator%iter_level = 0 406 iterator%ijk_nl = ijk_nl 407 408 CALL timestop(handle) 409 END SUBROUTINE 410 411! ************************************************************************************************** 412!> \brief Destroy 3c-nl iterator 413!> \param iterator ... 414! ************************************************************************************************** 415 SUBROUTINE neighbor_list_3c_iterator_destroy(iterator) 416 TYPE(neighbor_list_3c_iterator_type), & 417 INTENT(INOUT) :: iterator 418 419 CHARACTER(len=*), PARAMETER :: routineN = 'neighbor_list_3c_iterator_destroy', & 420 routineP = moduleN//':'//routineN 421 422 INTEGER :: handle 423 424 CALL timeset(routineN, handle) 425 CALL neighbor_list_iterator_release(iterator%iter_ij) 426 CALL neighbor_list_iterator_release(iterator%iter_jk) 427 NULLIFY (iterator%iter_ij) 428 NULLIFY (iterator%iter_jk) 429 430 CALL timestop(handle) 431 END SUBROUTINE 432 433! ************************************************************************************************** 434!> \brief Iterate 3c-nl iterator 435!> \param iterator ... 436!> \return 0 if successful; 1 if end was reached 437! ************************************************************************************************** 438 RECURSIVE FUNCTION neighbor_list_3c_iterate(iterator) RESULT(iter_stat) 439 TYPE(neighbor_list_3c_iterator_type), & 440 INTENT(INOUT) :: iterator 441 INTEGER :: iter_stat 442 443 CHARACTER(len=*), PARAMETER :: routineN = 'neighbor_list_3c_iterate', & 444 routineP = moduleN//':'//routineN 445 446 INTEGER :: iatom, iter_level, jatom, jatom_1, & 447 jatom_2, katom 448 LOGICAL :: skip_this 449 450 iter_level = iterator%iter_level 451 452 IF (iter_level == 0) THEN 453 iter_stat = neighbor_list_iterate(iterator%iter_ij) 454 455 IF (iter_stat /= 0) THEN 456 RETURN 457 ENDIF 458 ENDIF 459 iter_stat = nl_sub_iterate(iterator%iter_jk, iterator%iter_ij) 460 IF (iter_stat /= 0) THEN 461 iterator%iter_level = 0 462 iter_stat = neighbor_list_3c_iterate(iterator) 463 RETURN 464 ELSE 465 iterator%iter_level = 1 466 ENDIF 467 468 CPASSERT(iter_stat == 0) 469 CPASSERT(iterator%iter_level == 1) 470 CALL get_iterator_info(iterator%iter_ij, iatom=iatom, jatom=jatom_1) 471 CALL get_iterator_info(iterator%iter_jk, iatom=jatom_2, jatom=katom) 472 473 CPASSERT(jatom_1 == jatom_2) 474 jatom = jatom_1 475 476 skip_this = .TRUE. 477 SELECT CASE (iterator%ijk_nl%sym) 478 CASE (symmetric_none) 479 skip_this = .FALSE. 480 CASE (symmetric_ij) 481 skip_this = .NOT. include_symmetric(iatom, jatom) 482 CASE (symmetric_jk) 483 skip_this = .NOT. include_symmetric(jatom, katom) 484 CASE (symmetrik_ik) 485 skip_this = .NOT. include_symmetric(iatom, katom) 486 CASE (symmetric_ijk) 487 skip_this = .NOT. include_symmetric(iatom, jatom) .OR. .NOT. include_symmetric(jatom, katom) 488 CASE DEFAULT 489 CPABORT("should not happen") 490 END SELECT 491 492 IF (skip_this) THEN 493 iter_stat = neighbor_list_3c_iterate(iterator) 494 RETURN 495 ENDIF 496 497 END FUNCTION 498 499! ************************************************************************************************** 500!> \brief Get info of current iteration 501!> \param iterator ... 502!> \param ikind ... 503!> \param jkind ... 504!> \param kkind ... 505!> \param nkind ... 506!> \param iatom ... 507!> \param jatom ... 508!> \param katom ... 509!> \param rij ... 510!> \param rjk ... 511!> \param rik ... 512!> \param cell_j ... 513!> \param cell_k ... 514!> \return ... 515! ************************************************************************************************** 516 SUBROUTINE get_3c_iterator_info(iterator, ikind, jkind, kkind, nkind, iatom, jatom, katom, & 517 rij, rjk, rik, cell_j, cell_k) 518 TYPE(neighbor_list_3c_iterator_type), & 519 INTENT(INOUT) :: iterator 520 INTEGER, INTENT(OUT), OPTIONAL :: ikind, jkind, kkind, nkind, iatom, & 521 jatom, katom 522 REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: rij, rjk, rik 523 INTEGER, DIMENSION(3), INTENT(OUT), OPTIONAL :: cell_j, cell_k 524 525 CHARACTER(len=*), PARAMETER :: routineN = 'get_3c_iterator_info', & 526 routineP = moduleN//':'//routineN 527 528 INTEGER, DIMENSION(2) :: atoms_1, atoms_2, kinds_1, kinds_2 529 INTEGER, DIMENSION(3) :: cell_1, cell_2 530 REAL(KIND=dp), DIMENSION(3) :: r_1, r_2 531 532 CPASSERT(iterator%iter_level == 1) 533 534 CALL get_iterator_info(iterator%iter_ij, & 535 ikind=kinds_1(1), jkind=kinds_1(2), nkind=nkind, & 536 iatom=atoms_1(1), jatom=atoms_1(2), r=r_1, & 537 cell=cell_1) 538 539 CALL get_iterator_info(iterator%iter_jk, & 540 ikind=kinds_2(1), jkind=kinds_2(2), & 541 iatom=atoms_2(1), jatom=atoms_2(2), r=r_2, & 542 cell=cell_2) 543 544 IF (PRESENT(ikind)) ikind = kinds_1(1) 545 IF (PRESENT(jkind)) jkind = kinds_1(2) 546 IF (PRESENT(kkind)) kkind = kinds_2(2) 547 IF (PRESENT(iatom)) iatom = atoms_1(1) 548 IF (PRESENT(jatom)) jatom = atoms_1(2) 549 IF (PRESENT(katom)) katom = atoms_2(2) 550 551 IF (PRESENT(rij)) rij = r_1 552 IF (PRESENT(rjk)) rjk = r_2 553 IF (PRESENT(rik)) rik = r_1 + r_2 554 555 IF (PRESENT(cell_j)) cell_j = cell_1 556 IF (PRESENT(cell_k)) cell_k = cell_1 + cell_2 557 558 END SUBROUTINE 559 560! ************************************************************************************************** 561!> \brief Allocate blocks of a 3-center tensor based on neighborlist 562!> \param t3c empty DBCSR tensor 563!> Should be of shape (1,1) if no kpoints are used and of shape (nimages, nimages) 564!> if k-points are used 565!> \param nl_3c 3-center neighborlist 566!> \param basis_i ... 567!> \param basis_j ... 568!> \param basis_k ... 569!> \param qs_env ... 570!> \param potential_parameter ... 571!> \param op_pos ... 572!> \param do_kpoints ... 573! ************************************************************************************************** 574 SUBROUTINE alloc_block_3c(t3c, nl_3c, basis_i, basis_j, basis_k, qs_env, potential_parameter, op_pos, do_kpoints) 575 TYPE(dbcsr_t_type), DIMENSION(:, :), INTENT(INOUT) :: t3c 576 TYPE(neighbor_list_3c_type), INTENT(INOUT) :: nl_3c 577 TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j, basis_k 578 TYPE(qs_environment_type), POINTER :: qs_env 579 TYPE(libint_potential_type), INTENT(IN) :: potential_parameter 580 INTEGER, INTENT(IN), OPTIONAL :: op_pos 581 LOGICAL, INTENT(IN), OPTIONAL :: do_kpoints 582 583 CHARACTER(LEN=*), PARAMETER :: routineN = 'alloc_block_3c', routineP = moduleN//':'//routineN 584 585 INTEGER :: blk_cnt, handle, i, i_img, iatom, iblk, ikind, iproc, j_img, jatom, jcell, jkind, & 586 katom, kcell, kkind, natom, nimg, op_ij, op_jk, op_pos_prv 587 INTEGER, ALLOCATABLE, DIMENSION(:) :: tmp 588 INTEGER, DIMENSION(3) :: cell_j, cell_k 589 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 590 LOGICAL :: do_kpoints_prv, new_block 591 REAL(KIND=dp) :: dij, dik, djk, dr_ij, dr_ik, dr_jk, & 592 kind_radius_i, kind_radius_j, & 593 kind_radius_k 594 REAL(KIND=dp), DIMENSION(3) :: rij, rik, rjk 595 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 596 TYPE(cp_para_env_type), POINTER :: para_env 597 TYPE(dft_control_type), POINTER :: dft_control 598 TYPE(kpoint_type), POINTER :: kpoints 599 TYPE(neighbor_list_3c_iterator_type) :: nl_3c_iter 600 TYPE(one_dim_int_array), ALLOCATABLE, & 601 DIMENSION(:, :) :: alloc_i, alloc_j, alloc_k 602 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 603 604 CALL timeset(routineN, handle) 605 NULLIFY (qs_kind_set, atomic_kind_set) 606 607 IF (PRESENT(do_kpoints)) THEN 608 do_kpoints_prv = do_kpoints 609 ELSE 610 do_kpoints_prv = .FALSE. 611 ENDIF 612 613 dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp 614 615 op_ij = do_potential_id; op_jk = do_potential_id 616 617 IF (PRESENT(op_pos)) THEN 618 op_pos_prv = op_pos 619 ELSE 620 op_pos_prv = 1 621 ENDIF 622 623 SELECT CASE (op_pos_prv) 624 CASE (1) 625 op_ij = potential_parameter%potential_type 626 CASE (2) 627 op_jk = potential_parameter%potential_type 628 END SELECT 629 630 IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short) THEN 631 dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor 632 dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor 633 ELSEIF (op_ij == do_potential_coulomb) THEN 634 dr_ij = 1000000.0_dp 635 dr_ik = 1000000.0_dp 636 ENDIF 637 638 IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short) THEN 639 dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor 640 dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor 641 ELSEIF (op_jk == do_potential_coulomb) THEN 642 dr_jk = 1000000.0_dp 643 dr_ik = 1000000.0_dp 644 ENDIF 645 646 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, natom=natom, & 647 dft_control=dft_control, kpoints=kpoints, para_env=para_env) 648 649 IF (do_kpoints_prv) THEN 650 nimg = dft_control%nimages 651 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index) 652 ELSE 653 nimg = 1 654 END IF 655 656 ALLOCATE (alloc_i(nimg, nimg)) 657 ALLOCATE (alloc_j(nimg, nimg)) 658 ALLOCATE (alloc_k(nimg, nimg)) 659 660 CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c) 661 DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0) 662 CALL get_3c_iterator_info(nl_3c_iter, ikind=ikind, jkind=jkind, kkind=kkind, & 663 iatom=iatom, jatom=jatom, katom=katom, & 664 rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k) 665 666 IF (do_kpoints_prv) THEN 667 668 jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3)) 669 IF (jcell > nimg) CYCLE 670 671 kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3)) 672 IF (kcell > nimg) CYCLE 673 ELSE 674 jcell = 1; kcell = 1 675 END IF 676 677 djk = NORM2(rjk) 678 dij = NORM2(rij) 679 dik = NORM2(rik) 680 681 CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=kind_radius_i) 682 CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, kind_radius=kind_radius_j) 683 CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, kind_radius=kind_radius_k) 684 685 IF (kind_radius_j + kind_radius_i + dr_ij < dij) CYCLE 686 IF (kind_radius_j + kind_radius_k + dr_jk < djk) CYCLE 687 IF (kind_radius_k + kind_radius_i + dr_ik < dik) CYCLE 688 689 ! tensor is not symmetric therefore need to allocate rows and columns in 690 ! correspondence with neighborlist. Note that this only allocates half 691 ! of the blocks (since neighborlist is symmetric). After filling the blocks, 692 ! tensor will be added to its transposed 693 694 ASSOCIATE (ai=>alloc_i(jcell, kcell)) 695 ASSOCIATE (aj=>alloc_j(jcell, kcell)) 696 ASSOCIATE (ak=>alloc_k(jcell, kcell)) 697 698 new_block = .TRUE. 699 IF (ALLOCATED(aj%array)) THEN 700 DO iblk = 1, SIZE(aj%array) 701 IF (ai%array(iblk) == iatom .AND. & 702 aj%array(iblk) == jatom .AND. & 703 ak%array(iblk) == katom) THEN 704 new_block = .FALSE. 705 EXIT 706 ENDIF 707 ENDDO 708 ENDIF 709 IF (.NOT. new_block) CYCLE 710 711 IF (ALLOCATED(ai%array)) THEN 712 blk_cnt = SIZE(ai%array) 713 ALLOCATE (tmp(blk_cnt)) 714 tmp(:) = ai%array(:) 715 DEALLOCATE (ai%array) 716 ALLOCATE (ai%array(blk_cnt + 1)) 717 ai%array(1:blk_cnt) = tmp(:) 718 ai%array(blk_cnt + 1) = iatom 719 ELSE 720 ALLOCATE (ai%array(1)) 721 ai%array(1) = iatom 722 ENDIF 723 724 IF (ALLOCATED(aj%array)) THEN 725 tmp(:) = aj%array(:) 726 DEALLOCATE (aj%array) 727 ALLOCATE (aj%array(blk_cnt + 1)) 728 aj%array(1:blk_cnt) = tmp(:) 729 aj%array(blk_cnt + 1) = jatom 730 ELSE 731 ALLOCATE (aj%array(1)) 732 aj%array(1) = jatom 733 ENDIF 734 735 IF (ALLOCATED(ak%array)) THEN 736 tmp(:) = ak%array(:) 737 DEALLOCATE (ak%array) 738 ALLOCATE (ak%array(blk_cnt + 1)) 739 ak%array(1:blk_cnt) = tmp(:) 740 ak%array(blk_cnt + 1) = katom 741 ELSE 742 ALLOCATE (ak%array(1)) 743 ak%array(1) = katom 744 ENDIF 745 746 IF (ALLOCATED(tmp)) DEALLOCATE (tmp) 747 END ASSOCIATE 748 END ASSOCIATE 749 END ASSOCIATE 750 ENDDO 751 752 CALL neighbor_list_3c_iterator_destroy(nl_3c_iter) 753 754 DO i_img = 1, nimg 755 DO j_img = 1, nimg 756 IF (ALLOCATED(alloc_i(i_img, j_img)%array)) THEN 757 DO i = 1, SIZE(alloc_i(i_img, j_img)%array) 758 CALL dbcsr_t_get_stored_coordinates(t3c(i_img, j_img), & 759 [alloc_i(i_img, j_img)%array(i), alloc_j(i_img, j_img)%array(i), & 760 alloc_k(i_img, j_img)%array(i)], & 761 iproc) 762 CPASSERT(iproc .EQ. para_env%mepos) 763 ENDDO 764 765 CALL dbcsr_t_reserve_blocks(t3c(i_img, j_img), & 766 alloc_i(i_img, j_img)%array, & 767 alloc_j(i_img, j_img)%array, & 768 alloc_k(i_img, j_img)%array) 769 ENDIF 770 ENDDO 771 ENDDO 772 773 CALL timestop(handle) 774 775 END SUBROUTINE 776 777! ************************************************************************************************** 778!> \brief Build 3-center integral tensor 779!> \param t3c empty DBCSR tensor 780!> Should be of shape (1,1) if no kpoints are used and of shape (nimages, nimages) 781!> if k-points are used 782!> \param filter_eps Filter threshold for tensor blocks 783!> \param qs_env ... 784!> \param nl_3c 3-center neighborlist 785!> \param basis_i ... 786!> \param basis_j ... 787!> \param basis_k ... 788!> \param potential_parameter ... 789!> \param int_eps neglect integrals smaller than int_eps 790!> \param op_pos operator position. 791!> 1: calculate (i|jk) integrals, 792!> 2: calculate (ij|k) integrals 793!> \param do_kpoints ... 794!> this routine requires that libint has been static initialised somewhere else 795!> \param desymmetrize ... 796! ************************************************************************************************** 797 SUBROUTINE build_3c_integrals(t3c, filter_eps, qs_env, & 798 nl_3c, basis_i, basis_j, basis_k, & 799 potential_parameter, & 800 int_eps, & 801 op_pos, do_kpoints, desymmetrize) 802 TYPE(dbcsr_t_type), DIMENSION(:, :), INTENT(INOUT) :: t3c 803 REAL(KIND=dp), INTENT(IN) :: filter_eps 804 TYPE(qs_environment_type), POINTER :: qs_env 805 TYPE(neighbor_list_3c_type), INTENT(INOUT) :: nl_3c 806 TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j, basis_k 807 TYPE(libint_potential_type), INTENT(IN) :: potential_parameter 808 REAL(KIND=dp), INTENT(IN), OPTIONAL :: int_eps 809 INTEGER, INTENT(IN), OPTIONAL :: op_pos 810 LOGICAL, INTENT(IN), OPTIONAL :: do_kpoints, desymmetrize 811 812 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_3c_integrals', & 813 routineP = moduleN//':'//routineN 814 815 INTEGER :: block_end_i, block_end_j, block_end_k, block_start_i, block_start_j, & 816 block_start_k, handle, handle2, i, iatom, ibasis, ikind, imax, iset, jatom, jcell, jkind, & 817 jset, katom, kcell, kkind, kset, m_max, maxli, maxlj, maxlk, natom, ncoi, ncoj, ncok, & 818 nimg, nseti, nsetj, nsetk, op_ij, op_jk, op_pos_prv, sgfi, sgfj, sgfk, unit_id 819 INTEGER, DIMENSION(3) :: blk_size, cell_j, cell_k, sp 820 INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, & 821 lmin_k, npgfi, npgfj, npgfk, nsgfi, & 822 nsgfj, nsgfk 823 INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j, first_sgf_k 824 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 825 LOGICAL :: block_not_zero, debug, desymmetrize_prv, & 826 do_kpoints_prv, found 827 REAL(KIND=dp) :: dij, dik, djk, dr_ij, dr_ik, dr_jk, & 828 kind_radius_i, kind_radius_j, & 829 kind_radius_k, max_contraction_val, & 830 prefac 831 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: max_contraction_i, max_contraction_j, & 832 max_contraction_k 833 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: block_t, dummy_block_t, sijk, sijk_contr 834 REAL(KIND=dp), DIMENSION(3) :: ri, rij, rik, rj, rjk, rk 835 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j, set_radius_k 836 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, rpgf_k, sphi_i, sphi_j, & 837 sphi_k, zeti, zetj, zetk 838 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 839 TYPE(cp_libint_t) :: lib 840 TYPE(cp_para_env_type), POINTER :: para_env 841 TYPE(dbcsr_t_type) :: t_3c_tmp 842 TYPE(dft_control_type), POINTER :: dft_control 843 TYPE(kpoint_type), POINTER :: kpoints 844 TYPE(neighbor_list_3c_iterator_type) :: nl_3c_iter 845 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 846 847 CALL timeset(routineN, handle) 848 849 debug = .FALSE. 850 851 IF (PRESENT(do_kpoints)) THEN 852 do_kpoints_prv = do_kpoints 853 ELSE 854 do_kpoints_prv = .FALSE. 855 ENDIF 856 857 IF (PRESENT(desymmetrize)) THEN 858 desymmetrize_prv = desymmetrize 859 ELSE 860 desymmetrize_prv = .TRUE. 861 ENDIF 862 863 op_ij = do_potential_id; op_jk = do_potential_id 864 865 IF (PRESENT(op_pos)) THEN 866 op_pos_prv = op_pos 867 ELSE 868 op_pos_prv = 1 869 ENDIF 870 871 SELECT CASE (op_pos_prv) 872 CASE (1) 873 op_ij = potential_parameter%potential_type 874 CASE (2) 875 op_jk = potential_parameter%potential_type 876 END SELECT 877 878 dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp 879 880 IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short) THEN 881 dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor 882 dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor 883 ELSEIF (op_ij == do_potential_coulomb) THEN 884 dr_ij = 1000000.0_dp 885 dr_ik = 1000000.0_dp 886 ENDIF 887 888 IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short) THEN 889 dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor 890 dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor 891 ELSEIF (op_jk == do_potential_coulomb) THEN 892 dr_jk = 1000000.0_dp 893 dr_ik = 1000000.0_dp 894 ENDIF 895 896 NULLIFY (qs_kind_set, atomic_kind_set) 897 898 CALL alloc_block_3c(t3c, nl_3c, basis_i, basis_j, basis_k, qs_env, & 899 potential_parameter, op_pos=op_pos_prv, do_kpoints=do_kpoints) 900 901 ! get stuff 902 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, & 903 natom=natom, kpoints=kpoints, dft_control=dft_control, para_env=para_env) 904 905 IF (do_kpoints_prv) THEN 906 nimg = dft_control%nimages 907 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index) 908 ELSE 909 nimg = 1 910 END IF 911 912 CPASSERT(ALL(SHAPE(t3c) == [nimg, nimg])) 913 914 !Need the max l for each basis for libint 915 maxli = 0 916 DO ibasis = 1, SIZE(basis_i) 917 CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax) 918 maxli = MAX(maxli, imax) 919 END DO 920 maxlj = 0 921 DO ibasis = 1, SIZE(basis_j) 922 CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax) 923 maxlj = MAX(maxlj, imax) 924 END DO 925 maxlk = 0 926 DO ibasis = 1, SIZE(basis_k) 927 CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax) 928 maxlk = MAX(maxlk, imax) 929 END DO 930 m_max = maxli + maxlj + maxlk 931 932 !Init the truncated Coulomb operator 933 IF (op_ij == do_potential_truncated .OR. op_jk == do_potential_truncated) THEN 934 935 IF (m_max > get_lmax_init()) THEN 936 IF (para_env%mepos == 0) THEN 937 CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename) 938 END IF 939 CALL init(m_max, unit_id, para_env%mepos, para_env%group) 940 IF (para_env%mepos == 0) THEN 941 CALL close_file(unit_id) 942 END IF 943 END IF 944 END IF 945 946 CALL init_md_ftable(nmax=m_max) 947 948 IF (op_ij /= do_potential_id .OR. op_jk /= do_potential_id) THEN 949 CALL cp_libint_init_3eri(lib, MAX(maxli, maxlj, maxlk)) 950 CALL cp_libint_set_contrdepth(lib, 1) 951 ENDIF 952 953 CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c) 954 DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0) 955 CALL get_3c_iterator_info(nl_3c_iter, ikind=ikind, jkind=jkind, kkind=kkind, & 956 iatom=iatom, jatom=jatom, katom=katom, & 957 rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k) 958 959 IF (do_kpoints_prv) THEN 960 prefac = 0.5_dp 961 ELSEIF (nl_3c%sym == symmetric_jk) THEN 962 IF (jatom == katom) THEN 963 ! factor 0.5 due to double-counting of diagonal blocks 964 ! (we desymmetrize by adding transpose) 965 prefac = 0.5_dp 966 ELSE 967 prefac = 1.0_dp 968 ENDIF 969 ELSE 970 prefac = 1.0_dp 971 ENDIF 972 973 IF (do_kpoints_prv) THEN 974 975 jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3)) 976 IF (jcell > nimg) CYCLE 977 978 kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3)) 979 IF (kcell > nimg) CYCLE 980 981 ELSE 982 jcell = 1; kcell = 1 983 END IF 984 985 CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, & 986 npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, & 987 sphi=sphi_i, zet=zeti, kind_radius=kind_radius_i) 988 989 CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, & 990 npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, & 991 sphi=sphi_j, zet=zetj, kind_radius=kind_radius_j) 992 993 CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, first_sgf=first_sgf_k, lmax=lmax_k, lmin=lmin_k, & 994 npgf=npgfk, nset=nsetk, nsgf_set=nsgfk, pgf_radius=rpgf_k, set_radius=set_radius_k, & 995 sphi=sphi_k, zet=zetk, kind_radius=kind_radius_k) 996 997 djk = NORM2(rjk) 998 dij = NORM2(rij) 999 dik = NORM2(rik) 1000 1001 IF (kind_radius_j + kind_radius_i + dr_ij < dij) CYCLE 1002 IF (kind_radius_j + kind_radius_k + dr_jk < djk) CYCLE 1003 IF (kind_radius_k + kind_radius_i + dr_ik < dik) CYCLE 1004 1005 ALLOCATE (max_contraction_i(nseti)) 1006 max_contraction_i = 0.0_dp 1007 DO iset = 1, nseti 1008 sgfi = first_sgf_i(1, iset) 1009 max_contraction_i(iset) = MAXVAL((/(SUM(ABS(sphi_i(:, i))), i=sgfi, sgfi + nsgfi(iset) - 1)/)) 1010 ENDDO 1011 1012 ALLOCATE (max_contraction_j(nsetj)) 1013 max_contraction_j = 0.0_dp 1014 DO jset = 1, nsetj 1015 sgfj = first_sgf_j(1, jset) 1016 max_contraction_j(jset) = MAXVAL((/(SUM(ABS(sphi_j(:, i))), i=sgfj, sgfj + nsgfj(jset) - 1)/)) 1017 ENDDO 1018 1019 ALLOCATE (max_contraction_k(nsetk)) 1020 max_contraction_k = 0.0_dp 1021 DO kset = 1, nsetk 1022 sgfk = first_sgf_k(1, kset) 1023 max_contraction_k(kset) = MAXVAL((/(SUM(ABS(sphi_k(:, i))), i=sgfk, sgfk + nsgfk(kset) - 1)/)) 1024 ENDDO 1025 1026 CALL dbcsr_t_blk_sizes(t3c(jcell, kcell), [iatom, jatom, katom], blk_size) 1027 1028 IF (op_ij /= do_potential_id) THEN 1029 ALLOCATE (block_t(blk_size(2), blk_size(3), blk_size(1))) 1030 ELSE 1031 ALLOCATE (block_t(blk_size(1), blk_size(2), blk_size(3))) 1032 ENDIF 1033 1034 block_t = 0.0_dp 1035 block_not_zero = .FALSE. 1036 1037 DO iset = 1, nseti 1038 1039 DO jset = 1, nsetj 1040 1041 IF (set_radius_j(jset) + set_radius_i(iset) + dr_ij < dij) CYCLE 1042 1043 DO kset = 1, nsetk 1044 1045 IF (set_radius_j(jset) + set_radius_k(kset) + dr_jk < djk) CYCLE 1046 IF (set_radius_k(kset) + set_radius_i(iset) + dr_ik < dik) CYCLE 1047 1048 ncoi = npgfi(iset)*ncoset(lmax_i(iset)) 1049 ncoj = npgfj(jset)*ncoset(lmax_j(jset)) 1050 ncok = npgfk(kset)*ncoset(lmax_k(kset)) 1051 1052 sgfi = first_sgf_i(1, iset) 1053 sgfj = first_sgf_j(1, jset) 1054 sgfk = first_sgf_k(1, kset) 1055 1056 IF (ncoj*ncok*ncoi > 0) THEN 1057 1058 IF (op_ij == do_potential_id .AND. op_jk == do_potential_id) THEN ! cp2k overlap integrals 1059 ALLOCATE (sijk(ncoi, ncoj, ncok)) 1060 sijk(:, :, :) = 0.0_dp 1061 CALL overlap3(lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), & 1062 lmin_i(iset), & 1063 lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), lmin_j(jset), & 1064 lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), lmin_k(kset), & 1065 rij, dij, rik, dik, rjk, djk, sijk) 1066 ELSEIF (op_jk /= do_potential_id) THEN ! for everything else use libint 1067 ALLOCATE (sijk(ncoi, ncoj, ncok)) 1068 sijk(:, :, :) = 0.0_dp 1069 !need positions for libint. Only relative positions are needed => set ri to 0.0 1070 ri = 0.0_dp 1071 rj = rij ! ri + rij 1072 rk = rik ! ri + rik 1073 CALL eri_3center(sijk, & 1074 lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), ri, & 1075 lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), rj, & 1076 lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), rk, & 1077 dij, dik, djk, lib, potential_parameter) 1078 ELSEIF (op_ij /= do_potential_id) THEN 1079 ALLOCATE (sijk(ncoj, ncok, ncoi)) 1080 sijk(:, :, :) = 0.0_dp 1081 ri = 0.0_dp 1082 rj = rij ! ri + rij 1083 rk = rik ! ri + rik 1084 CALL eri_3center(sijk, & 1085 lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), rj, & 1086 lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), rk, & 1087 lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), ri, & 1088 djk, dij, dik, lib, potential_parameter) 1089 1090 ENDIF 1091 max_contraction_val = max_contraction_i(iset)*max_contraction_j(jset)*max_contraction_k(kset) 1092 1093 IF (PRESENT(int_eps)) THEN 1094 IF (MAXVAL(ABS(sijk))*max_contraction_val < int_eps) THEN 1095 DEALLOCATE (sijk) 1096 CYCLE 1097 ENDIF 1098 ENDIF 1099 1100 block_not_zero = .TRUE. 1101 1102 IF (op_ij /= do_potential_id) THEN 1103 ALLOCATE (sijk_contr(nsgfj(jset), nsgfk(kset), nsgfi(iset))) 1104 CALL abc_contract(sijk_contr, sijk, & 1105 sphi_j(:, sgfj:), sphi_k(:, sgfk:), sphi_i(:, sgfi:), & 1106 ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), nsgfi(iset)) 1107 DEALLOCATE (sijk) 1108 1109 block_start_j = sgfj 1110 block_end_j = sgfj + nsgfj(jset) - 1 1111 block_start_k = sgfk 1112 block_end_k = sgfk + nsgfk(kset) - 1 1113 block_start_i = sgfi 1114 block_end_i = sgfi + nsgfi(iset) - 1 1115 1116 block_t(block_start_j:block_end_j, & 1117 block_start_k:block_end_k, & 1118 block_start_i:block_end_i) = & 1119 block_t(block_start_j:block_end_j, & 1120 block_start_k:block_end_k, & 1121 block_start_i:block_end_i) + & 1122 prefac*sijk_contr(:, :, :) 1123 DEALLOCATE (sijk_contr) 1124 ELSE 1125 ALLOCATE (sijk_contr(nsgfi(iset), nsgfj(jset), nsgfk(kset))) 1126 1127 CALL abc_contract(sijk_contr, sijk, & 1128 sphi_i(:, sgfi:), sphi_j(:, sgfj:), sphi_k(:, sgfk:), & 1129 ncoi, ncoj, ncok, nsgfi(iset), nsgfj(jset), nsgfk(kset)) 1130 1131 DEALLOCATE (sijk) 1132 block_start_j = sgfj 1133 block_end_j = sgfj + nsgfj(jset) - 1 1134 block_start_k = sgfk 1135 block_end_k = sgfk + nsgfk(kset) - 1 1136 block_start_i = sgfi 1137 block_end_i = sgfi + nsgfi(iset) - 1 1138 1139 block_t(block_start_i:block_end_i, & 1140 block_start_j:block_end_j, & 1141 block_start_k:block_end_k) = & 1142 block_t(block_start_i:block_end_i, & 1143 block_start_j:block_end_j, & 1144 block_start_k:block_end_k) + & 1145 prefac*sijk_contr(:, :, :) 1146 DEALLOCATE (sijk_contr) 1147 1148 ENDIF 1149 1150 END IF ! number of triples > 0 1151 1152 END DO 1153 1154 END DO 1155 1156 END DO 1157 1158 IF (block_not_zero) THEN 1159 CALL timeset(routineN//"_put_dbcsr", handle2) 1160 IF (debug) THEN 1161 CALL dbcsr_t_get_block(t3c(jcell, kcell), & 1162 [iatom, jatom, katom], dummy_block_t, found=found) 1163 CPASSERT(found) 1164 ENDIF 1165 sp = SHAPE(block_t) 1166 1167 IF (op_ij /= do_potential_id) THEN 1168 1169 sp([2, 3, 1]) = sp 1170 1171 CALL dbcsr_t_put_block(t3c(jcell, kcell), & 1172 [iatom, jatom, katom], sp, RESHAPE(block_t, SHAPE=sp, ORDER=[2, 3, 1]), summation=.TRUE.) 1173 ELSE 1174 CALL dbcsr_t_put_block(t3c(jcell, kcell), & 1175 [iatom, jatom, katom], sp, block_t, summation=.TRUE.) 1176 ENDIF 1177 CALL timestop(handle2) 1178 ENDIF 1179 1180 DEALLOCATE (block_t) 1181 1182 DEALLOCATE (max_contraction_i, max_contraction_j, max_contraction_k) 1183 END DO 1184 1185 IF (op_ij /= do_potential_id .OR. op_jk /= do_potential_id) THEN 1186 CALL cp_libint_cleanup_3eri(lib) 1187 ENDIF 1188 1189 CALL neighbor_list_3c_iterator_destroy(nl_3c_iter) 1190 1191 IF (nl_3c%sym == symmetric_jk .OR. do_kpoints_prv) THEN 1192 DO jcell = 1, nimg 1193 DO kcell = 1, nimg 1194 ! need half of filter eps because afterwards we add transposed tensor 1195 CALL dbcsr_t_filter(t3c(jcell, kcell), filter_eps/2) 1196 ENDDO 1197 ENDDO 1198 1199 IF (desymmetrize_prv) THEN 1200 ! add transposed of overlap integrals 1201 CALL dbcsr_t_create(t3c(1, 1), t_3c_tmp) 1202 DO jcell = 1, nimg 1203 DO kcell = 1, jcell 1204 CALL dbcsr_t_copy(t3c(jcell, kcell), t_3c_tmp) 1205 CALL dbcsr_t_copy(t_3c_tmp, t3c(kcell, jcell), order=[1, 3, 2], summation=.TRUE., move_data=.TRUE.) 1206 CALL dbcsr_t_filter(t3c(kcell, jcell), filter_eps) 1207 ENDDO 1208 ENDDO 1209 DO jcell = 1, nimg 1210 DO kcell = jcell + 1, nimg 1211 CALL dbcsr_t_copy(t3c(jcell, kcell), t_3c_tmp) 1212 CALL dbcsr_t_copy(t_3c_tmp, t3c(kcell, jcell), order=[1, 3, 2], summation=.FALSE., move_data=.TRUE.) 1213 CALL dbcsr_t_filter(t3c(kcell, jcell), filter_eps) 1214 ENDDO 1215 ENDDO 1216 CALL dbcsr_t_destroy(t_3c_tmp) 1217 ENDIF 1218 ELSEIF (nl_3c%sym == symmetric_none) THEN 1219 DO jcell = 1, nimg 1220 DO kcell = 1, nimg 1221 CALL dbcsr_t_filter(t3c(jcell, kcell), filter_eps) 1222 ENDDO 1223 ENDDO 1224 ELSE 1225 CPABORT("requested symmetric case not implemented") 1226 ENDIF 1227 1228 CALL timestop(handle) 1229 END SUBROUTINE 1230 1231! ************************************************************************************************** 1232!> \brief contiguous distribution of weighted elements 1233!> \param nel ... 1234!> \param nbin ... 1235!> \param weights ... 1236!> \param limits_start ... 1237!> \param limits_end ... 1238!> \param dist ... 1239! ************************************************************************************************** 1240 SUBROUTINE contiguous_tensor_dist(nel, nbin, weights, limits_start, limits_end, dist) 1241 INTEGER, INTENT(IN) :: nel, nbin 1242 INTEGER, DIMENSION(nel), INTENT(IN) :: weights 1243 INTEGER, DIMENSION(nbin), INTENT(OUT), OPTIONAL :: limits_start, limits_end 1244 INTEGER, DIMENSION(nel), INTENT(OUT), OPTIONAL :: dist 1245 1246 INTEGER :: el_end, el_start, end_weight, ibin, & 1247 nel_div, nel_rem, nel_split, nel_w, & 1248 w_partialsum 1249 1250 nel_w = SUM(weights) 1251 nel_div = nel_w/nbin 1252 nel_rem = MOD(nel_w, nbin) 1253 1254 w_partialsum = 0 1255 el_end = 0 1256 end_weight = 0 1257 DO ibin = 1, nbin 1258 nel_split = nel_div 1259 IF (ibin <= nel_rem) THEN 1260 nel_split = nel_split + 1 1261 ENDIF 1262 el_start = el_end + 1 1263 el_end = el_start 1264 w_partialsum = w_partialsum + weights(el_end) 1265 end_weight = end_weight + nel_split 1266 DO WHILE (w_partialsum < end_weight) 1267 !IF (ABS(w_partialsum + weights(el_end) - end_weight) > ABS(w_partialsum - end_weight)) EXIT 1268 el_end = el_end + 1 1269 w_partialsum = w_partialsum + weights(el_end) 1270 ENDDO 1271 IF (PRESENT(dist)) dist(el_start:el_end) = ibin - 1 1272 IF (PRESENT(limits_start)) limits_start(ibin) = el_start 1273 IF (PRESENT(limits_end)) limits_end(ibin) = el_end 1274 ENDDO 1275 1276 END SUBROUTINE contiguous_tensor_dist 1277 1278! ************************************************************************************************** 1279!> \brief ... 1280!> \param t2c empty DBCSR matrix 1281!> \param filter_eps Filter threshold for matrix blocks 1282!> \param qs_env ... 1283!> \param nl_2c 2-center neighborlist 1284!> \param basis_i ... 1285!> \param basis_j ... 1286!> \param potential_parameter ... 1287!> \param do_kpoints ... 1288!> this routine requires that libint has been static initialised somewhere else 1289! ************************************************************************************************** 1290 SUBROUTINE build_2c_integrals(t2c, filter_eps, qs_env, & 1291 nl_2c, basis_i, basis_j, & 1292 potential_parameter, do_kpoints) 1293 TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: t2c 1294 REAL(KIND=dp), INTENT(IN) :: filter_eps 1295 TYPE(qs_environment_type), POINTER :: qs_env 1296 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 1297 POINTER :: nl_2c 1298 TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j 1299 TYPE(libint_potential_type), INTENT(IN) :: potential_parameter 1300 LOGICAL, INTENT(IN), OPTIONAL :: do_kpoints 1301 1302 CHARACTER(len=*), PARAMETER :: routineN = 'build_2c_integrals', & 1303 routineP = moduleN//':'//routineN 1304 1305 INTEGER :: handle, iatom, ibasis, icol, ikind, imax, img, irow, iset, jatom, jkind, jset, & 1306 m_max, maxli, maxlj, n1, n2, natom, ncoi, ncoj, nimg, nseti, nsetj, offi, offj, op_prv, & 1307 sgfi, sgfj, unit_id 1308 INTEGER, DIMENSION(3) :: cell 1309 INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmin_i, lmin_j, npgfi, & 1310 npgfj, nsgfi, nsgfj 1311 INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j 1312 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 1313 LOGICAL :: do_kpoints_prv, do_symmetric, found, & 1314 trans 1315 REAL(KIND=dp) :: dab 1316 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sij, sij_contr, sij_rs 1317 REAL(KIND=dp), DIMENSION(3) :: ri, rij, rj 1318 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j 1319 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, scon_i, scon_j, sphi_i, & 1320 sphi_j, zeti, zetj 1321 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1322 TYPE(block_p_type) :: block_t 1323 TYPE(cp_libint_t) :: lib 1324 TYPE(cp_para_env_type), POINTER :: para_env 1325 TYPE(dft_control_type), POINTER :: dft_control 1326 TYPE(kpoint_type), POINTER :: kpoints 1327 TYPE(neighbor_list_iterator_p_type), & 1328 DIMENSION(:), POINTER :: nl_iterator 1329 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1330 1331 CALL timeset(routineN, handle) 1332 1333 IF (PRESENT(do_kpoints)) THEN 1334 do_kpoints_prv = do_kpoints 1335 ELSE 1336 do_kpoints_prv = .FALSE. 1337 ENDIF 1338 1339 op_prv = potential_parameter%potential_type 1340 1341 NULLIFY (qs_kind_set, atomic_kind_set, block_t%block, cell_to_index) 1342 1343 ! get stuff 1344 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, & 1345 natom=natom, kpoints=kpoints, dft_control=dft_control, para_env=para_env) 1346 1347 IF (do_kpoints_prv) THEN 1348 nimg = dft_control%nimages 1349 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index) 1350 ELSE 1351 nimg = 1 1352 END IF 1353 1354 CPASSERT(ALL(SHAPE(t2c) == [nimg])) 1355 1356 ! check for symmetry 1357 CPASSERT(SIZE(nl_2c) > 0) 1358 CALL get_neighbor_list_set_p(neighbor_list_sets=nl_2c, symmetric=do_symmetric) 1359 1360 IF (do_symmetric) THEN 1361 DO img = 1, nimg 1362 CPASSERT(dbcsr_has_symmetry(t2c(img))) 1363 ENDDO 1364 ELSE 1365 DO img = 1, nimg 1366 CPASSERT(.NOT. dbcsr_has_symmetry(t2c(img))) 1367 ENDDO 1368 ENDIF 1369 1370 DO img = 1, nimg 1371 CALL cp_dbcsr_alloc_block_from_nbl(t2c(img), nl_2c) 1372 ENDDO 1373 1374 maxli = 0 1375 DO ibasis = 1, SIZE(basis_i) 1376 CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax) 1377 maxli = MAX(maxli, imax) 1378 END DO 1379 maxlj = 0 1380 DO ibasis = 1, SIZE(basis_j) 1381 CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax) 1382 maxlj = MAX(maxlj, imax) 1383 END DO 1384 1385 m_max = maxli + maxlj 1386 1387 !Init the truncated Coulomb operator 1388 IF (op_prv == do_potential_truncated) THEN 1389 1390 IF (m_max > get_lmax_init()) THEN 1391 IF (para_env%mepos == 0) THEN 1392 CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename) 1393 END IF 1394 CALL init(m_max, unit_id, para_env%mepos, para_env%group) 1395 IF (para_env%mepos == 0) THEN 1396 CALL close_file(unit_id) 1397 END IF 1398 END IF 1399 END IF 1400 1401 CALL init_md_ftable(nmax=m_max) 1402 1403 IF (op_prv /= do_potential_id) THEN 1404 CALL cp_libint_init_2eri(lib, MAX(maxli, maxlj)) 1405 CALL cp_libint_set_contrdepth(lib, 1) 1406 ENDIF 1407 1408 CALL neighbor_list_iterator_create(nl_iterator, nl_2c) 1409 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 1410 1411 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, & 1412 iatom=iatom, jatom=jatom, r=rij, cell=cell) 1413 IF (do_kpoints_prv) THEN 1414 img = cell_to_index(cell(1), cell(2), cell(3)) 1415 IF (img > nimg) CYCLE 1416 ELSE 1417 img = 1 1418 END IF 1419 1420 CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, & 1421 npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, & 1422 sphi=sphi_i, zet=zeti, scon=scon_i) 1423 1424 CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, & 1425 npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, & 1426 sphi=sphi_j, zet=zetj, scon=scon_j) 1427 1428 IF (do_symmetric) THEN 1429 IF (iatom <= jatom) THEN 1430 irow = iatom 1431 icol = jatom 1432 ELSE 1433 irow = jatom 1434 icol = iatom 1435 END IF 1436 ELSE 1437 irow = iatom 1438 icol = jatom 1439 END IF 1440 1441 dab = NORM2(rij) 1442 1443 CALL dbcsr_get_block_p(matrix=t2c(img), & 1444 row=irow, col=icol, BLOCK=block_t%block, found=found) 1445 CPASSERT(found) 1446 trans = do_symmetric .AND. (iatom > jatom) 1447 1448 DO iset = 1, nseti 1449 1450 ncoi = npgfi(iset)*ncoset(lmax_i(iset)) 1451 n1 = npgfi(iset)*(ncoset(lmax_i(iset)) - ncoset(lmin_i(iset) - 1)) 1452 sgfi = first_sgf_i(1, iset) 1453 offi = ncoset(lmin_i(iset) - 1) + 1 1454 1455 DO jset = 1, nsetj 1456 1457 ncoj = npgfj(jset)*ncoset(lmax_j(jset)) 1458 n2 = npgfj(jset)*(ncoset(lmax_j(jset)) - ncoset(lmin_j(jset) - 1)) 1459 sgfj = first_sgf_j(1, jset) 1460 offj = ncoset(lmin_j(jset) - 1) + 1 1461 1462 IF (ncoi*ncoj > 0) THEN 1463 ALLOCATE (sij_contr(nsgfi(iset), nsgfj(jset))) 1464 sij_contr(:, :) = 0.0_dp 1465 1466 IF (op_prv == do_potential_id) THEN 1467 ALLOCATE (sij(n1, n2)) 1468 sij(:, :) = 0.0_dp 1469 1470 CALL overlap_ab(lmax_i(iset), lmin_i(iset), npgfi(iset), rpgf_i(:, iset), zeti(:, iset), & 1471 lmax_j(jset), lmin_j(jset), npgfj(jset), rpgf_j(:, jset), zetj(:, jset), & 1472 rij, sab=sij(:, :)) 1473 1474 CALL ab_contract(sij_contr, sij, & 1475 scon_i(:, sgfi:), scon_j(:, sgfj:), & 1476 n1, n2, nsgfi(iset), nsgfj(jset)) 1477 1478 ELSE 1479 ALLOCATE (sij(ncoi, ncoj)) 1480 sij(:, :) = 0.0_dp 1481 1482 ri = 0.0_dp 1483 rj = rij 1484 1485 CALL eri_2center(sij, lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), & 1486 rpgf_i(:, iset), ri, lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), & 1487 rpgf_j(:, jset), rj, dab, lib, potential_parameter) 1488 1489 CALL ab_contract(sij_contr, sij, & 1490 sphi_i(:, sgfi:), sphi_j(:, sgfj:), & 1491 ncoi, ncoj, nsgfi(iset), nsgfj(jset)) 1492 ENDIF 1493 1494 DEALLOCATE (sij) 1495 IF (trans) THEN 1496 ALLOCATE (sij_rs(nsgfj(jset), nsgfi(iset))) 1497 sij_rs(:, :) = TRANSPOSE(sij_contr) 1498 ELSE 1499 ALLOCATE (sij_rs(nsgfi(iset), nsgfj(jset))) 1500 sij_rs(:, :) = sij_contr 1501 ENDIF 1502 1503 DEALLOCATE (sij_contr) 1504 1505 CALL block_add("IN", sij_rs, & 1506 nsgfi(iset), nsgfj(jset), block_t%block, & 1507 sgfi, sgfj, trans=trans) 1508 DEALLOCATE (sij_rs) 1509 ENDIF 1510 END DO 1511 END DO 1512 ENDDO 1513 1514 IF (op_prv /= do_potential_id) THEN 1515 CALL cp_libint_cleanup_2eri(lib) 1516 ENDIF 1517 1518 CALL neighbor_list_iterator_release(nl_iterator) 1519 DO img = 1, nimg 1520 CALL dbcsr_finalize(t2c(img)) 1521 CALL dbcsr_filter(t2c(img), filter_eps) 1522 ENDDO 1523 1524 CALL timestop(handle) 1525 1526 END SUBROUTINE 1527 1528! ************************************************************************************************** 1529!> \brief Change process grid of tensor, a load balanced default distribution adapted to the new grid 1530!> is chosen automatically. 1531!> \param t_3c ... 1532!> \param pgrid ... 1533!> \param nodata ... 1534!> \param unit_nr ... 1535! ************************************************************************************************** 1536 SUBROUTINE tensor_change_pgrid(t_3c, pgrid, nodata, unit_nr) 1537 TYPE(dbcsr_t_type), INTENT(INOUT) :: t_3c 1538 TYPE(dbcsr_t_pgrid_type), INTENT(IN) :: pgrid 1539 LOGICAL, INTENT(IN), OPTIONAL :: nodata 1540 INTEGER, INTENT(IN), OPTIONAL :: unit_nr 1541 1542 CHARACTER(LEN=*), PARAMETER :: routineN = 'tensor_change_pgrid', & 1543 routineP = moduleN//':'//routineN 1544 1545 CHARACTER(default_string_length) :: name 1546 INTEGER :: handle 1547 INTEGER, ALLOCATABLE, DIMENSION(:) :: bs1, bs2, bs3, dist1, dist2, dist3, & 1548 map1, map2 1549 INTEGER, DIMENSION(3) :: pcoord, pcoord_ref, pdims, pdims_ref, & 1550 tdims 1551 TYPE(dbcsr_t_distribution_type) :: dist 1552 TYPE(dbcsr_t_type) :: t_tmp 1553 1554 CALL dbcsr_t_mp_environ_pgrid(pgrid, pdims, pcoord) 1555 CALL dbcsr_t_mp_environ_pgrid(t_3c%pgrid, pdims_ref, pcoord_ref) 1556 1557 IF (ALL(pdims == pdims_ref)) RETURN 1558 1559 CALL timeset(routineN, handle) 1560 1561 CALL dbcsr_t_get_info(t_3c, nblks_total=tdims, blk_size_1=bs1, blk_size_2=bs2, blk_size_3=bs3, & 1562 name=name) 1563 1564 ALLOCATE (dist1(tdims(1))) 1565 ALLOCATE (dist2(tdims(2))) 1566 ALLOCATE (dist3(tdims(3))) 1567 1568 CALL cyclic_tensor_dist(tdims(1), pdims(1), bs1, dist1) 1569 CALL cyclic_tensor_dist(tdims(2), pdims(2), bs2, dist2) 1570 CALL cyclic_tensor_dist(tdims(3), pdims(3), bs3, dist3) 1571 1572 CALL dbcsr_t_get_mapping_info(t_3c%nd_index_blk, map1_2d=map1, map2_2d=map2) 1573 CALL dbcsr_t_distribution_new(dist, pgrid, map1, map2, dist1, dist2, dist3) 1574 CALL dbcsr_t_create(t_tmp, name, dist, map1, map2, dbcsr_type_real_8, bs1, bs2, bs3) 1575 CALL dbcsr_t_distribution_destroy(dist) 1576 1577 IF (PRESENT(nodata)) THEN 1578 IF (.NOT. nodata) CALL dbcsr_t_copy(t_3c, t_tmp, move_data=.TRUE.) 1579 ELSE 1580 CALL dbcsr_t_copy(t_3c, t_tmp, move_data=.TRUE.) 1581 ENDIF 1582 1583 CALL dbcsr_t_destroy(t_3c) 1584 t_3c = t_tmp 1585 1586 IF (PRESENT(unit_nr)) THEN 1587 IF (unit_nr > 0) THEN 1588 WRITE (unit_nr, "(T2,A,1X,A)") "OPTIMIZED PGRID INFO FOR", TRIM(t_3c%name) 1589 WRITE (unit_nr, "(T4,A,1X,3I6)") "process grid dimensions:", pdims 1590 CALL dbcsr_t_write_split_info(pgrid, unit_nr) 1591 ENDIF 1592 ENDIF 1593 1594 CALL timestop(handle) 1595 END SUBROUTINE 1596 1597END MODULE 1598