1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Generate the atomic neighbor lists. 8!> \par History 9!> - List rebuild for sab_orb neighbor list (10.09.2002,MK) 10!> - List rebuild for all lists (25.09.2002,MK) 11!> - Row-wise parallelized version (16.06.2003,MK) 12!> - Row- and column-wise parallelized version (19.07.2003,MK) 13!> - bug fix for non-periodic case (23.02.06,MK) 14!> - major refactoring (25.07.10,jhu) 15!> \author Matthias Krack (08.10.1999,26.03.2002,16.06.2003) 16! ************************************************************************************************** 17MODULE qs_neighbor_lists 18 USE almo_scf_types, ONLY: almo_max_cutoff_multiplier 19 USE atomic_kind_types, ONLY: atomic_kind_type,& 20 get_atomic_kind,& 21 get_atomic_kind_set 22 USE basis_set_types, ONLY: get_gto_basis_set,& 23 gto_basis_set_p_type,& 24 gto_basis_set_type 25 USE cell_types, ONLY: cell_type,& 26 get_cell,& 27 pbc,& 28 plane_distance,& 29 real_to_scaled,& 30 scaled_to_real 31 USE cp_control_types, ONLY: dft_control_type 32 USE cp_log_handling, ONLY: cp_get_default_logger,& 33 cp_logger_type 34 USE cp_output_handling, ONLY: cp_p_file,& 35 cp_print_key_finished_output,& 36 cp_print_key_should_output,& 37 cp_print_key_unit_nr 38 USE cp_para_types, ONLY: cp_para_env_type 39 USE cp_units, ONLY: cp_unit_from_cp2k 40 USE distribution_1d_types, ONLY: distribution_1d_type 41 USE distribution_2d_types, ONLY: distribution_2d_type 42 USE ewald_environment_types, ONLY: ewald_env_get,& 43 ewald_environment_type 44 USE external_potential_types, ONLY: all_potential_type,& 45 get_potential,& 46 gth_potential_type,& 47 sgp_potential_type 48 USE input_constants, ONLY: dispersion_uff,& 49 do_method_lrigpw,& 50 do_method_rigpw,& 51 do_se_IS_slater,& 52 vdw_pairpot_dftd3,& 53 vdw_pairpot_dftd3bj,& 54 xc_vdw_fun_pairpot 55 USE input_section_types, ONLY: section_vals_get,& 56 section_vals_get_subs_vals,& 57 section_vals_type,& 58 section_vals_val_get 59 USE kinds, ONLY: default_string_length,& 60 dp,& 61 int_8 62 USE kpoint_types, ONLY: kpoint_type 63 USE message_passing, ONLY: mp_max,& 64 mp_sum 65 USE molecule_types, ONLY: molecule_type 66 USE particle_types, ONLY: particle_type 67 USE paw_proj_set_types, ONLY: get_paw_proj_set,& 68 paw_proj_set_type 69 USE periodic_table, ONLY: ptable 70 USE physcon, ONLY: bohr 71 USE qs_dftb_types, ONLY: qs_dftb_atom_type 72 USE qs_dftb_utils, ONLY: get_dftb_atom_param 73 USE qs_dispersion_types, ONLY: qs_dispersion_type 74 USE qs_environment_types, ONLY: get_qs_env,& 75 qs_environment_type 76 USE qs_gcp_types, ONLY: qs_gcp_type 77 USE qs_kind_types, ONLY: get_qs_kind,& 78 get_qs_kind_set,& 79 qs_kind_type 80 USE qs_ks_types, ONLY: get_ks_env,& 81 qs_ks_env_type,& 82 set_ks_env 83 USE qs_neighbor_list_types, ONLY: & 84 add_neighbor_list, add_neighbor_node, allocate_neighbor_list_set, & 85 deallocate_neighbor_list_set, get_iterator_info, neighbor_list_iterate, & 86 neighbor_list_iterator_create, neighbor_list_iterator_p_type, & 87 neighbor_list_iterator_release, neighbor_list_p_type, neighbor_list_set_p_type, & 88 neighbor_list_set_type 89 USE string_utilities, ONLY: compress 90 USE subcell_types, ONLY: allocate_subcell,& 91 deallocate_subcell,& 92 give_ijk_subcell,& 93 subcell_type 94 USE util, ONLY: locate,& 95 sort 96#include "./base/base_uses.f90" 97 98 IMPLICIT NONE 99 100 PRIVATE 101 102! ************************************************************************************************** 103 TYPE local_atoms_type 104 INTEGER, DIMENSION(:), POINTER :: list, & 105 list_local_a_index, & 106 list_local_b_index, & 107 list_1d, & 108 list_a_mol, & 109 list_b_mol 110 END TYPE local_atoms_type 111! ************************************************************************************************** 112 113 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_neighbor_lists' 114 115 ! private counter, used to version qs neighbor lists 116 INTEGER, SAVE, PRIVATE :: last_qs_neighbor_list_id_nr = 0 117 118 ! Public subroutines 119 PUBLIC :: build_qs_neighbor_lists, local_atoms_type, atom2d_cleanup, & 120 atom2d_build, build_neighbor_lists, pair_radius_setup, & 121 setup_neighbor_list, write_neighbor_lists 122CONTAINS 123 124! ************************************************************************************************** 125!> \brief free the internals of atom2d 126!> \param atom2d ... 127!> \param 128! ************************************************************************************************** 129 SUBROUTINE atom2d_cleanup(atom2d) 130 TYPE(local_atoms_type), DIMENSION(:) :: atom2d 131 132 CHARACTER(len=*), PARAMETER :: routineN = 'atom2d_cleanup', routineP = moduleN//':'//routineN 133 134 INTEGER :: handle, ikind 135 136 CALL timeset(routineN, handle) 137 DO ikind = 1, SIZE(atom2d) 138 NULLIFY (atom2d(ikind)%list) 139 IF (ASSOCIATED(atom2d(ikind)%list_local_a_index)) THEN 140 DEALLOCATE (atom2d(ikind)%list_local_a_index) 141 END IF 142 IF (ASSOCIATED(atom2d(ikind)%list_local_b_index)) THEN 143 DEALLOCATE (atom2d(ikind)%list_local_b_index) 144 END IF 145 IF (ASSOCIATED(atom2d(ikind)%list_a_mol)) THEN 146 DEALLOCATE (atom2d(ikind)%list_a_mol) 147 END IF 148 IF (ASSOCIATED(atom2d(ikind)%list_b_mol)) THEN 149 DEALLOCATE (atom2d(ikind)%list_b_mol) 150 END IF 151 IF (ASSOCIATED(atom2d(ikind)%list_1d)) THEN 152 DEALLOCATE (atom2d(ikind)%list_1d) 153 END IF 154 END DO 155 CALL timestop(handle) 156 157 END SUBROUTINE atom2d_cleanup 158 159! ************************************************************************************************** 160!> \brief Build some distribution structure of atoms, refactored from build_qs_neighbor_lists 161!> \param atom2d output 162!> \param distribution_1d ... 163!> \param distribution_2d ... 164!> \param atomic_kind_set ... 165!> \param molecule_set ... 166!> \param molecule_only ... 167!> \param particle_set ... 168!> \author JH 169! ************************************************************************************************** 170 SUBROUTINE atom2d_build(atom2d, distribution_1d, distribution_2d, & 171 atomic_kind_set, molecule_set, molecule_only, particle_set) 172 TYPE(local_atoms_type), DIMENSION(:) :: atom2d 173 TYPE(distribution_1d_type), POINTER :: distribution_1d 174 TYPE(distribution_2d_type), POINTER :: distribution_2d 175 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 176 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 177 LOGICAL :: molecule_only 178 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 179 180 CHARACTER(len=*), PARAMETER :: routineN = 'atom2d_build', routineP = moduleN//':'//routineN 181 182 INTEGER :: atom_a, handle, ia, iat, iatom, & 183 iatom_local, ikind, imol, natom, & 184 natom_a, natom_local_a, natom_local_b, & 185 nel, nkind 186 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom2mol, listindex, listsort 187 INTEGER, DIMENSION(:), POINTER :: atom_of_kind, local_cols_array, & 188 local_rows_array 189 190 CALL timeset(routineN, handle) 191 192 nkind = SIZE(atomic_kind_set) 193 natom = SIZE(particle_set) 194 ALLOCATE (atom_of_kind(natom)) 195 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 196 atom_of_kind=atom_of_kind) 197 198 IF (molecule_only) THEN 199 ALLOCATE (atom2mol(natom)) 200 DO imol = 1, SIZE(molecule_set) 201 DO iat = molecule_set(imol)%first_atom, molecule_set(imol)%last_atom 202 atom2mol(iat) = imol 203 ENDDO 204 END DO 205 ENDIF 206 207 DO ikind = 1, nkind 208 NULLIFY (atom2d(ikind)%list) 209 NULLIFY (atom2d(ikind)%list_local_a_index) 210 NULLIFY (atom2d(ikind)%list_local_b_index) 211 NULLIFY (atom2d(ikind)%list_1d) 212 NULLIFY (atom2d(ikind)%list_a_mol) 213 NULLIFY (atom2d(ikind)%list_b_mol) 214 215 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom2d(ikind)%list) 216 217 natom_a = SIZE(atom2d(ikind)%list) 218 219 natom_local_a = distribution_2d%n_local_rows(ikind) 220 natom_local_b = distribution_2d%n_local_cols(ikind) 221 local_rows_array => distribution_2d%local_rows(ikind)%array 222 local_cols_array => distribution_2d%local_cols(ikind)%array 223 224 nel = distribution_1d%n_el(ikind) 225 ALLOCATE (atom2d(ikind)%list_1d(nel)) 226 DO iat = 1, nel 227 ia = distribution_1d%list(ikind)%array(iat) 228 atom2d(ikind)%list_1d(iat) = atom_of_kind(ia) 229 END DO 230 231 ALLOCATE (listsort(natom_a), listindex(natom_a)) 232 listsort(1:natom_a) = atom2d(ikind)%list(1:natom_a) 233 CALL sort(listsort, natom_a, listindex) 234 ! Block rows 235 IF (natom_local_a > 0) THEN 236 ALLOCATE (atom2d(ikind)%list_local_a_index(natom_local_a)) 237 ALLOCATE (atom2d(ikind)%list_a_mol(natom_local_a)) 238 atom2d(ikind)%list_a_mol(:) = 0 239 240 ! Build index vector for mapping 241 DO iatom_local = 1, natom_local_a 242 atom_a = local_rows_array(iatom_local) 243 iatom = locate(listsort, atom_a) 244 atom2d(ikind)%list_local_a_index(iatom_local) = listindex(iatom) 245 IF (molecule_only) atom2d(ikind)%list_a_mol(iatom_local) = atom2mol(atom_a) 246 END DO 247 248 END IF 249 250 ! Block columns 251 IF (natom_local_b > 0) THEN 252 253 ALLOCATE (atom2d(ikind)%list_local_b_index(natom_local_b)) 254 ALLOCATE (atom2d(ikind)%list_b_mol(natom_local_b)) 255 atom2d(ikind)%list_b_mol(:) = 0 256 257 ! Build index vector for mapping 258 DO iatom_local = 1, natom_local_b 259 atom_a = local_cols_array(iatom_local) 260 iatom = locate(listsort, atom_a) 261 atom2d(ikind)%list_local_b_index(iatom_local) = listindex(iatom) 262 IF (molecule_only) atom2d(ikind)%list_b_mol(iatom_local) = atom2mol(atom_a) 263 END DO 264 265 END IF 266 267 DEALLOCATE (listsort, listindex) 268 269 ENDDO 270 271 DEALLOCATE (atom_of_kind) 272 273 CALL timestop(handle) 274 275 END SUBROUTINE atom2d_build 276 277! ************************************************************************************************** 278!> \brief Build all the required neighbor lists for Quickstep. 279!> \param qs_env ... 280!> \param para_env ... 281!> \param molecular ... 282!> \param force_env_section ... 283!> \date 28.08.2000 284!> \par History 285!> - Major refactoring (25.07.2010,jhu) 286!> \author MK 287!> \version 1.0 288! ************************************************************************************************** 289 SUBROUTINE build_qs_neighbor_lists(qs_env, para_env, molecular, force_env_section) 290 TYPE(qs_environment_type), POINTER :: qs_env 291 TYPE(cp_para_env_type), POINTER :: para_env 292 LOGICAL, OPTIONAL :: molecular 293 TYPE(section_vals_type), POINTER :: force_env_section 294 295 CHARACTER(len=*), PARAMETER :: routineN = 'build_qs_neighbor_lists', & 296 routineP = moduleN//':'//routineN 297 298 CHARACTER(LEN=default_string_length) :: print_key_path 299 INTEGER :: handle, ikind, iw, maxatom, nkind, & 300 numshells, zat 301 LOGICAL :: all_potential_present, almo, dftb, do_hfx, dokp, explicit, gth_potential_present, & 302 lri_optbas, lrigpw, mic, molecule_only, nddo, paw_atom, paw_atom_present, rigpw, & 303 sgp_potential_present, xtb 304 LOGICAL, ALLOCATABLE, DIMENSION(:) :: all_present, aux_fit_present, aux_present, & 305 core_present, default_present, oce_present, orb_present, ppl_present, ppnl_present, & 306 ri_present, xb1_atom, xb2_atom 307 REAL(dp) :: almo_rcov, almo_rvdw, rcut, roperator, & 308 subcells 309 REAL(dp), ALLOCATABLE, DIMENSION(:) :: all_pot_rad, aux_fit_radius, c_radius, calpha, & 310 core_radius, oce_radius, orb_radius, ppl_radius, ppnl_radius, ri_radius, zeff 311 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius 312 TYPE(all_potential_type), POINTER :: all_potential 313 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 314 TYPE(cell_type), POINTER :: cell 315 TYPE(cp_logger_type), POINTER :: logger 316 TYPE(dft_control_type), POINTER :: dft_control 317 TYPE(distribution_1d_type), POINTER :: distribution_1d 318 TYPE(distribution_2d_type), POINTER :: distribution_2d 319 TYPE(ewald_environment_type), POINTER :: ewald_env 320 TYPE(gth_potential_type), POINTER :: gth_potential 321 TYPE(gto_basis_set_type), POINTER :: aux_basis_set, aux_fit_basis_set, & 322 orb_basis_set, ri_basis_set 323 TYPE(kpoint_type), POINTER :: kpoints 324 TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d 325 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 326 TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: saa_list, sab_all, sab_almo, & 327 sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, sab_cn, sab_core, sab_gcp, sab_kp, & 328 sab_lrc, sab_orb, sab_scp, sab_se, sab_tbe, sab_vdw, sab_xb, sac_ae, sac_lri, sac_ppl, & 329 sap_oce, sap_ppnl, soa_list, soo_list 330 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 331 TYPE(paw_proj_set_type), POINTER :: paw_proj 332 TYPE(qs_dftb_atom_type), POINTER :: dftb_atom 333 TYPE(qs_dispersion_type), POINTER :: dispersion_env 334 TYPE(qs_gcp_type), POINTER :: gcp_env 335 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 336 TYPE(qs_ks_env_type), POINTER :: ks_env 337 TYPE(section_vals_type), POINTER :: hfx_sections, neighbor_list_section 338 TYPE(sgp_potential_type), POINTER :: sgp_potential 339 340 CALL timeset(routineN, handle) 341 NULLIFY (logger) 342 logger => cp_get_default_logger() 343 344 NULLIFY (atomic_kind_set, qs_kind_set, cell, neighbor_list_section, & 345 distribution_1d, distribution_2d, gth_potential, sgp_potential, orb_basis_set, & 346 particle_set, molecule_set, dft_control, ks_env) 347 348 NULLIFY (sab_orb) 349 NULLIFY (sac_ae) 350 NULLIFY (sac_ppl) 351 NULLIFY (sac_lri) 352 NULLIFY (sap_ppnl) 353 NULLIFY (sap_oce) 354 NULLIFY (sab_se) 355 NULLIFY (sab_lrc) 356 NULLIFY (sab_tbe) 357 NULLIFY (sab_core) 358 NULLIFY (sab_xb) 359 NULLIFY (sab_all) 360 NULLIFY (sab_vdw) 361 NULLIFY (sab_cn) 362 NULLIFY (sab_aux_fit) 363 NULLIFY (sab_aux_fit_vs_orb) 364 NULLIFY (soo_list) 365 NULLIFY (sab_scp) 366 NULLIFY (sab_almo) 367 NULLIFY (sab_kp) 368 369 CALL get_qs_env(qs_env, & 370 ks_env=ks_env, & 371 atomic_kind_set=atomic_kind_set, & 372 qs_kind_set=qs_kind_set, & 373 cell=cell, & 374 kpoints=kpoints, & 375 distribution_2d=distribution_2d, & 376 local_particles=distribution_1d, & 377 particle_set=particle_set, & 378 molecule_set=molecule_set, & 379 dft_control=dft_control) 380 381 neighbor_list_section => section_vals_get_subs_vals(force_env_section, "DFT%PRINT%NEIGHBOR_LISTS") 382 383 ! This sets the id number of the qs neighbor lists, new lists, means new version 384 ! new version implies new sparsity of the matrices 385 last_qs_neighbor_list_id_nr = last_qs_neighbor_list_id_nr + 1 386 CALL set_ks_env(ks_env=ks_env, neighbor_list_id=last_qs_neighbor_list_id_nr) 387 388 CALL get_ks_env(ks_env=ks_env, & 389 sab_orb=sab_orb, & 390 sab_aux_fit=sab_aux_fit, & 391 sab_aux_fit_vs_orb=sab_aux_fit_vs_orb, & 392 sab_aux_fit_asymm=sab_aux_fit_asymm, & 393 sac_ae=sac_ae, & 394 sac_ppl=sac_ppl, & 395 sac_lri=sac_lri, & 396 sab_vdw=sab_vdw, & 397 sap_ppnl=sap_ppnl, & 398 sap_oce=sap_oce, & 399 sab_se=sab_se, & 400 sab_lrc=sab_lrc, & 401 sab_tbe=sab_tbe, & 402 sab_core=sab_core, & 403 sab_xb=sab_xb, & 404 sab_scp=sab_scp, & 405 sab_all=sab_all, & 406 sab_almo=sab_almo, & 407 sab_kp=sab_kp) 408 409 dokp = (kpoints%nkp > 0) 410 nddo = dft_control%qs_control%semi_empirical 411 dftb = dft_control%qs_control%dftb 412 xtb = dft_control%qs_control%xtb 413 almo = dft_control%qs_control%do_almo_scf 414 lrigpw = (dft_control%qs_control%method_id == do_method_lrigpw) 415 rigpw = (dft_control%qs_control%method_id == do_method_rigpw) 416 lri_optbas = dft_control%qs_control%lri_optbas 417 418 ! molecular lists 419 molecule_only = .FALSE. 420 IF (PRESENT(molecular)) molecule_only = molecular 421 ! minimum image convention (MIC) 422 mic = molecule_only 423 IF (dokp) THEN 424 ! no MIC for kpoints 425 mic = .FALSE. 426 ELSEIF (nddo) THEN 427 ! enforce MIC for interaction lists in SE 428 mic = .TRUE. 429 END IF 430 431 hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF") 432 CALL section_vals_get(hfx_sections, explicit=do_hfx) 433 434 CALL get_atomic_kind_set(atomic_kind_set, maxatom=maxatom) 435 CALL get_qs_kind_set(qs_kind_set, paw_atom_present=paw_atom_present, & 436 gth_potential_present=gth_potential_present, & 437 sgp_potential_present=sgp_potential_present, & 438 all_potential_present=all_potential_present) 439 440 CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells) 441 442 ! Allocate work storage 443 nkind = SIZE(atomic_kind_set) 444 ALLOCATE (orb_present(nkind), aux_fit_present(nkind), aux_present(nkind), & 445 default_present(nkind), core_present(nkind)) 446 ALLOCATE (orb_radius(nkind), aux_fit_radius(nkind), c_radius(nkind), & 447 core_radius(nkind), calpha(nkind), zeff(nkind)) 448 orb_radius(:) = 0.0_dp 449 aux_fit_radius(:) = 0.0_dp 450 c_radius(:) = 0.0_dp 451 core_radius(:) = 0.0_dp 452 calpha(:) = 0.0_dp 453 zeff(:) = 0.0_dp 454 455 ALLOCATE (pair_radius(nkind, nkind)) 456 IF (gth_potential_present .OR. sgp_potential_present) THEN 457 ALLOCATE (ppl_present(nkind), ppl_radius(nkind)) 458 ppl_radius = 0.0_dp 459 ALLOCATE (ppnl_present(nkind), ppnl_radius(nkind)) 460 ppnl_radius = 0.0_dp 461 END IF 462 IF (paw_atom_present) THEN 463 ALLOCATE (oce_present(nkind), oce_radius(nkind)) 464 oce_radius = 0.0_dp 465 END IF 466 IF (all_potential_present .OR. sgp_potential_present) THEN 467 ALLOCATE (all_present(nkind), all_pot_rad(nkind)) 468 all_pot_rad = 0.0_dp 469 END IF 470 471 ! Initialize the local data structures 472 ALLOCATE (atom2d(nkind)) 473 CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, & 474 molecule_set, molecule_only, particle_set=particle_set) 475 476 DO ikind = 1, nkind 477 478 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom2d(ikind)%list) 479 480 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB") 481 CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_basis_set, basis_type="AUX") 482 CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_fit_basis_set, basis_type="AUX_FIT") 483 484 CALL get_qs_kind(qs_kind_set(ikind), & 485 paw_proj_set=paw_proj, & 486 paw_atom=paw_atom, & 487 all_potential=all_potential, & 488 gth_potential=gth_potential, & 489 sgp_potential=sgp_potential) 490 491 IF (dftb) THEN 492 ! Set the interaction radius for the neighbor lists (DFTB case) 493 ! This includes all interactions (orbitals and short range pair potential) except vdW 494 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom) 495 CALL get_dftb_atom_param(dftb_parameter=dftb_atom, & 496 cutoff=orb_radius(ikind), & 497 defined=orb_present(ikind)) 498 ELSE 499 IF (ASSOCIATED(orb_basis_set)) THEN 500 orb_present(ikind) = .TRUE. 501 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, kind_radius=orb_radius(ikind)) 502 ELSE 503 orb_present(ikind) = .FALSE. 504 ENDIF 505 ENDIF 506 507 IF (ASSOCIATED(aux_basis_set)) THEN 508 aux_present(ikind) = .TRUE. 509 ELSE 510 aux_present(ikind) = .FALSE. 511 ENDIF 512 513 IF (ASSOCIATED(aux_fit_basis_set)) THEN 514 aux_fit_present(ikind) = .TRUE. 515 CALL get_gto_basis_set(gto_basis_set=aux_fit_basis_set, kind_radius=aux_fit_radius(ikind)) 516 ELSE 517 aux_fit_present(ikind) = .FALSE. 518 END IF 519 520 ! core overlap 521 CALL get_qs_kind(qs_kind_set(ikind), & 522 alpha_core_charge=calpha(ikind), & 523 core_charge_radius=core_radius(ikind), & 524 zeff=zeff(ikind)) 525 IF (zeff(ikind) /= 0._dp .AND. calpha(ikind) /= 0._dp) THEN 526 core_present(ikind) = .TRUE. 527 ELSE 528 core_present(ikind) = .FALSE. 529 END IF 530 531 ! Pseudopotentials 532 IF (gth_potential_present .OR. sgp_potential_present) THEN 533 IF (ASSOCIATED(gth_potential)) THEN 534 CALL get_potential(potential=gth_potential, & 535 ppl_present=ppl_present(ikind), & 536 ppl_radius=ppl_radius(ikind), & 537 ppnl_present=ppnl_present(ikind), & 538 ppnl_radius=ppnl_radius(ikind)) 539 ELSE IF (ASSOCIATED(sgp_potential)) THEN 540 CALL get_potential(potential=sgp_potential, & 541 ppl_present=ppl_present(ikind), & 542 ppl_radius=ppl_radius(ikind), & 543 ppnl_present=ppnl_present(ikind), & 544 ppnl_radius=ppnl_radius(ikind)) 545 ELSE 546 ppl_present(ikind) = .FALSE. 547 ppnl_present(ikind) = .FALSE. 548 END IF 549 END IF 550 551 ! GAPW 552 IF (paw_atom_present) THEN 553 IF (paw_atom) THEN 554 oce_present(ikind) = .TRUE. 555 CALL get_paw_proj_set(paw_proj_set=paw_proj, rcprj=oce_radius(ikind)) 556 ELSE 557 oce_present(ikind) = .FALSE. 558 END IF 559 END IF 560 561 ! Check the presence of an all electron potential or ERFC potential 562 IF (all_potential_present .OR. sgp_potential_present) THEN 563 all_present(ikind) = .FALSE. 564 all_pot_rad(ikind) = 0.0_dp 565 IF (ASSOCIATED(all_potential)) THEN 566 all_present(ikind) = .TRUE. 567 CALL get_potential(potential=all_potential, core_charge_radius=all_pot_rad(ikind)) 568 ELSE IF (ASSOCIATED(sgp_potential)) THEN 569 IF (sgp_potential%ecp_local) THEN 570 all_present(ikind) = .TRUE. 571 CALL get_potential(potential=sgp_potential, core_charge_radius=all_pot_rad(ikind)) 572 END IF 573 END IF 574 END IF 575 576 END DO 577 578 ! Build the orbital-orbital overlap neighbor lists 579 CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius) 580 CALL build_neighbor_lists(sab_orb, particle_set, atom2d, cell, pair_radius, & 581 mic=mic, subcells=subcells, molecular=molecule_only, nlname="sab_orb") 582 CALL set_ks_env(ks_env=ks_env, sab_orb=sab_orb) 583 CALL write_neighbor_lists(sab_orb, particle_set, cell, para_env, neighbor_list_section, & 584 "/SAB_ORB", "sab_orb", "ORBITAL ORBITAL") 585 586 ! Build orbital-orbital list containing all the pairs, to be used with 587 ! non-symmetric operators. Beware: the cutoff of the orbital-orbital overlap 588 ! might not be optimal. It should be verified for each operator. 589 IF (.NOT. (nddo .OR. dftb .OR. xtb)) THEN 590 CALL build_neighbor_lists(sab_all, particle_set, atom2d, cell, pair_radius, & 591 mic=mic, symmetric=.FALSE., subcells=subcells, molecular=molecule_only, nlname="sab_all") 592 CALL set_ks_env(ks_env=ks_env, sab_all=sab_all) 593 ENDIF 594 595 ! Build the core-core overlap neighbor lists 596 IF (.NOT. (nddo .OR. dftb .OR. xtb)) THEN 597 CALL pair_radius_setup(core_present, core_present, core_radius, core_radius, pair_radius) 598 CALL build_neighbor_lists(sab_core, particle_set, atom2d, cell, pair_radius, subcells=subcells, & 599 operator_type="PP", nlname="sab_core") 600 CALL set_ks_env(ks_env=ks_env, sab_core=sab_core) 601 CALL write_neighbor_lists(sab_core, particle_set, cell, para_env, neighbor_list_section, & 602 "/SAB_CORE", "sab_core", "CORE CORE") 603 ENDIF 604 605 IF (dft_control%do_admm) THEN 606 CALL pair_radius_setup(aux_fit_present, aux_fit_present, aux_fit_radius, aux_fit_radius, pair_radius) 607 CALL build_neighbor_lists(sab_aux_fit, particle_set, atom2d, cell, pair_radius, & 608 mic=mic, molecular=molecule_only, subcells=subcells, nlname="sab_aux_fit") 609 CALL build_neighbor_lists(sab_aux_fit_asymm, particle_set, atom2d, cell, pair_radius, & 610 mic=mic, symmetric=.FALSE., molecular=molecule_only, subcells=subcells, & 611 nlname="sab_aux_fit_asymm") 612 CALL pair_radius_setup(aux_fit_present, orb_present, aux_fit_radius, orb_radius, pair_radius) 613 CALL build_neighbor_lists(sab_aux_fit_vs_orb, particle_set, atom2d, cell, pair_radius, & 614 mic=mic, symmetric=.FALSE., molecular=molecule_only, subcells=subcells, & 615 nlname="sab_aux_fit_vs_orb") 616 617 CALL set_ks_env(ks_env=ks_env, sab_aux_fit=sab_aux_fit) 618 CALL set_ks_env(ks_env=ks_env, sab_aux_fit_vs_orb=sab_aux_fit_vs_orb) 619 CALL set_ks_env(ks_env=ks_env, sab_aux_fit_asymm=sab_aux_fit_asymm) 620 621 CALL write_neighbor_lists(sab_aux_fit, particle_set, cell, para_env, neighbor_list_section, & 622 "/SAB_AUX_FIT", "sab_aux_fit", "AUX_FIT_ORBITAL AUX_FIT_ORBITAL") 623 CALL write_neighbor_lists(sab_aux_fit_vs_orb, particle_set, cell, para_env, neighbor_list_section, & 624 "/SAB_AUX_FIT_VS_ORB", "sab_aux_fit_vs_orb", "ORBITAL AUX_FIT_ORBITAL") 625 END IF 626 627 IF (dokp) THEN 628 ! We try to guess an integration radius for K-points 629 ! For non-HFX calculations we use the overlap list 630 ! For HFX we use the interaction radius of kinds (ORB or ADMM basis) 631 ! plus a range for the operator (some arbitrary guess, should be revisited!) 632 IF (do_hfx) THEN 633 roperator = 0.0_dp 634 CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%CUTOFF_RADIUS", & 635 explicit=explicit) 636 IF (explicit) THEN 637 CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%CUTOFF_RADIUS", & 638 r_val=roperator) 639 ELSE 640 CALL section_vals_val_get(hfx_sections, "PERIODIC%NUMBER_OF_SHELLS", & 641 i_val=numshells) 642 numshells = MAX(1, numshells) 643 roperator = MAX(plane_distance(1, 0, 0, cell), & 644 plane_distance(0, 1, 0, cell), & 645 plane_distance(0, 0, 1, cell))*numshells 646 END IF 647 IF (dft_control%do_admm) THEN 648 CALL pair_radius_setup(aux_fit_present, aux_fit_present, aux_fit_radius, aux_fit_radius, & 649 pair_radius) 650 ELSE 651 CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius) 652 END IF 653 pair_radius = pair_radius + roperator 654 ELSE 655 CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius) 656 END IF 657 CALL build_neighbor_lists(sab_kp, particle_set, atom2d, cell, pair_radius, & 658 subcells=subcells, nlname="sab_kp") 659 CALL set_ks_env(ks_env=ks_env, sab_kp=sab_kp) 660 END IF 661 662 ! Build orbital GTH-PPL operator overlap list 663 IF (gth_potential_present .OR. sgp_potential_present) THEN 664 IF (ANY(ppl_present)) THEN 665 CALL pair_radius_setup(orb_present, ppl_present, orb_radius, ppl_radius, pair_radius) 666 CALL build_neighbor_lists(sac_ppl, particle_set, atom2d, cell, pair_radius, & 667 subcells=subcells, operator_type="ABC", nlname="sac_ppl") 668 CALL set_ks_env(ks_env=ks_env, sac_ppl=sac_ppl) 669 CALL write_neighbor_lists(sac_ppl, particle_set, cell, para_env, neighbor_list_section, & 670 "/SAC_PPL", "sac_ppl", "ORBITAL GTH-PPL") 671 IF (lrigpw) THEN 672 IF (qs_env%lri_env%ppl_ri) THEN 673 CALL build_neighbor_lists(sac_lri, particle_set, atom2d, cell, pair_radius, & 674 subcells=subcells, symmetric=.FALSE., operator_type="PP", nlname="sac_lri") 675 CALL set_ks_env(ks_env=ks_env, sac_lri=sac_lri) 676 END IF 677 END IF 678 END IF 679 680 IF (ANY(ppnl_present)) THEN 681 CALL pair_radius_setup(orb_present, ppnl_present, orb_radius, ppnl_radius, pair_radius) 682 CALL build_neighbor_lists(sap_ppnl, particle_set, atom2d, cell, pair_radius, & 683 subcells=subcells, operator_type="ABBA", nlname="sap_ppnl") 684 CALL set_ks_env(ks_env=ks_env, sap_ppnl=sap_ppnl) 685 CALL write_neighbor_lists(sap_ppnl, particle_set, cell, para_env, neighbor_list_section, & 686 "/SAP_PPNL", "sap_ppnl", "ORBITAL GTH-PPNL") 687 END IF 688 END IF 689 690 IF (paw_atom_present) THEN 691 ! Build orbital-GAPW projector overlap list 692 IF (ANY(oce_present)) THEN 693 CALL pair_radius_setup(orb_present, oce_present, orb_radius, oce_radius, pair_radius) 694 CALL build_neighbor_lists(sap_oce, particle_set, atom2d, cell, pair_radius, & 695 subcells=subcells, operator_type="ABBA", nlname="sap_oce") 696 CALL set_ks_env(ks_env=ks_env, sap_oce=sap_oce) 697 CALL write_neighbor_lists(sap_oce, particle_set, cell, para_env, neighbor_list_section, & 698 "/SAP_OCE", "sap_oce", "ORBITAL(A) PAW-PRJ") 699 END IF 700 END IF 701 702 ! Build orbital-ERFC potential list 703 IF (.NOT. (nddo .OR. dftb .OR. xtb)) THEN 704 IF (all_potential_present .OR. sgp_potential_present) THEN 705 CALL pair_radius_setup(orb_present, all_present, orb_radius, all_pot_rad, pair_radius) 706 CALL build_neighbor_lists(sac_ae, particle_set, atom2d, cell, pair_radius, & 707 subcells=subcells, operator_type="ABC", nlname="sac_ae") 708 CALL set_ks_env(ks_env=ks_env, sac_ae=sac_ae) 709 CALL write_neighbor_lists(sac_ae, particle_set, cell, para_env, neighbor_list_section, & 710 "/SAC_AE", "sac_ae", "ORBITAL ERFC POTENTIAL") 711 END IF 712 END IF 713 714 IF (nddo) THEN 715 ! Semi-empirical neighbor lists 716 default_present = .TRUE. 717 c_radius = dft_control%qs_control%se_control%cutoff_cou 718 ! Build the neighbor lists for the Hartree terms 719 CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius) 720 IF (dft_control%qs_control%se_control%do_ewald_gks) THEN 721 ! Use MIC for the periodic code of GKS 722 CALL build_neighbor_lists(sab_se, particle_set, atom2d, cell, pair_radius, mic=mic, & 723 subcells=subcells, nlname="sab_se") 724 ELSE 725 CALL build_neighbor_lists(sab_se, particle_set, atom2d, cell, pair_radius, & 726 subcells=subcells, nlname="sab_se") 727 END IF 728 CALL set_ks_env(ks_env=ks_env, sab_se=sab_se) 729 CALL write_neighbor_lists(sab_se, particle_set, cell, para_env, neighbor_list_section, & 730 "/SAB_SE", "sab_se", "HARTREE INTERACTIONS") 731 732 ! If requested build the SE long-range correction neighbor list 733 IF ((dft_control%qs_control%se_control%do_ewald) .AND. & 734 (dft_control%qs_control%se_control%integral_screening /= do_se_IS_slater)) THEN 735 c_radius = dft_control%qs_control%se_control%cutoff_lrc 736 CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius) 737 CALL build_neighbor_lists(sab_lrc, particle_set, atom2d, cell, pair_radius, & 738 subcells=subcells, nlname="sab_lrc") 739 CALL set_ks_env(ks_env=ks_env, sab_lrc=sab_lrc) 740 CALL write_neighbor_lists(sab_lrc, particle_set, cell, para_env, neighbor_list_section, & 741 "/SAB_LRC", "sab_lrc", "SE LONG-RANGE CORRECTION") 742 END IF 743 END IF 744 745 IF (dftb) THEN 746 ! Build the neighbor lists for the DFTB Ewald methods 747 IF (dft_control%qs_control%dftb_control%do_ewald) THEN 748 CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env) 749 CALL ewald_env_get(ewald_env, rcut=rcut) 750 c_radius = rcut 751 CALL pair_radius_setup(orb_present, orb_present, c_radius, c_radius, pair_radius) 752 CALL build_neighbor_lists(sab_tbe, particle_set, atom2d, cell, pair_radius, mic=mic, & 753 subcells=subcells, nlname="sab_tbe") 754 CALL set_ks_env(ks_env=ks_env, sab_tbe=sab_tbe) 755 END IF 756 757 ! Build the neighbor lists for the DFTB vdW pair potential 758 IF (dft_control%qs_control%dftb_control%dispersion) THEN 759 IF (dft_control%qs_control%dftb_control%dispersion_type == dispersion_uff) THEN 760 DO ikind = 1, nkind 761 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom) 762 CALL get_dftb_atom_param(dftb_parameter=dftb_atom, rcdisp=c_radius(ikind)) 763 END DO 764 default_present = .TRUE. 765 CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius) 766 CALL build_neighbor_lists(sab_vdw, particle_set, atom2d, cell, pair_radius, & 767 subcells=subcells, nlname="sab_vdw") 768 CALL set_ks_env(ks_env=ks_env, sab_vdw=sab_vdw) 769 END IF 770 END IF 771 END IF 772 773 IF (xtb) THEN 774 ! Build the neighbor lists for the xTB Ewald method 775 IF (dft_control%qs_control%xtb_control%do_ewald) THEN 776 CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env) 777 CALL ewald_env_get(ewald_env, rcut=rcut) 778 c_radius = rcut 779 CALL pair_radius_setup(orb_present, orb_present, c_radius, c_radius, pair_radius) 780 CALL build_neighbor_lists(sab_tbe, particle_set, atom2d, cell, pair_radius, mic=mic, & 781 subcells=subcells, nlname="sab_tbe") 782 CALL set_ks_env(ks_env=ks_env, sab_tbe=sab_tbe) 783 END IF 784 ! XB list 785 ALLOCATE (xb1_atom(nkind), xb2_atom(nkind)) 786 c_radius = 0.5_dp*dft_control%qs_control%xtb_control%xb_radius 787 DO ikind = 1, nkind 788 CALL get_atomic_kind(atomic_kind_set(ikind), z=zat) 789 IF (zat == 17 .OR. zat == 35 .OR. zat == 53 .OR. zat == 85) THEN 790 xb1_atom(ikind) = .TRUE. 791 ELSE 792 xb1_atom(ikind) = .FALSE. 793 END IF 794 IF (zat == 7 .OR. zat == 8 .OR. zat == 15 .OR. zat == 16) THEN 795 xb2_atom(ikind) = .TRUE. 796 ELSE 797 xb2_atom(ikind) = .FALSE. 798 END IF 799 END DO 800 CALL pair_radius_setup(xb1_atom, xb2_atom, c_radius, c_radius, pair_radius) 801 CALL build_neighbor_lists(sab_xb, particle_set, atom2d, cell, pair_radius, & 802 symmetric=.FALSE., subcells=subcells, operator_type="PP", nlname="sab_xb") 803 CALL set_ks_env(ks_env=ks_env, sab_xb=sab_xb) 804 END IF 805 806 ! Build the neighbor lists for the vdW pair potential 807 CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env) 808 sab_vdw => dispersion_env%sab_vdw 809 sab_cn => dispersion_env%sab_cn 810 IF (dispersion_env%type == xc_vdw_fun_pairpot) THEN 811 c_radius(:) = dispersion_env%rc_disp 812 default_present = .TRUE. !include all atoms in vdW (even without basis) 813 CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius) 814 CALL build_neighbor_lists(sab_vdw, particle_set, atom2d, cell, pair_radius, & 815 subcells=subcells, operator_type="PP", nlname="sab_vdw") 816 dispersion_env%sab_vdw => sab_vdw 817 818 IF (dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. & 819 dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN 820 ! Build the neighbor lists for coordination numbers as needed by the DFT-D3 method 821 DO ikind = 1, nkind 822 CALL get_atomic_kind(atomic_kind_set(ikind), z=zat) 823 c_radius(ikind) = 4._dp*ptable(zat)%covalent_radius*bohr 824 END DO 825 CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius) 826 CALL build_neighbor_lists(sab_cn, particle_set, atom2d, cell, pair_radius, & 827 subcells=subcells, operator_type="PP", nlname="sab_cn") 828 dispersion_env%sab_cn => sab_cn 829 END IF 830 END IF 831 832 ! Build the neighbor lists for the gCP pair potential 833 NULLIFY (gcp_env) 834 CALL get_qs_env(qs_env=qs_env, gcp_env=gcp_env) 835 IF (ASSOCIATED(gcp_env)) THEN 836 IF (gcp_env%do_gcp) THEN 837 sab_gcp => gcp_env%sab_gcp 838 DO ikind = 1, nkind 839 c_radius(ikind) = gcp_env%gcp_kind(ikind)%rcsto 840 END DO 841 CALL pair_radius_setup(orb_present, orb_present, c_radius, c_radius, pair_radius) 842 CALL build_neighbor_lists(sab_gcp, particle_set, atom2d, cell, pair_radius, & 843 subcells=subcells, operator_type="PP", nlname="sab_gcp") 844 gcp_env%sab_gcp => sab_gcp 845 ELSE 846 NULLIFY (gcp_env%sab_gcp) 847 END IF 848 END IF 849 850 IF (lrigpw .OR. lri_optbas) THEN 851 ! set neighborlists in lri_env environment 852 CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius) 853 soo_list => qs_env%lri_env%soo_list 854 CALL build_neighbor_lists(soo_list, particle_set, atom2d, cell, pair_radius, & 855 mic=mic, molecular=molecule_only, subcells=subcells, nlname="soo_list") 856 qs_env%lri_env%soo_list => soo_list 857 CALL write_neighbor_lists(soo_list, particle_set, cell, para_env, neighbor_list_section, & 858 "/SOO_LIST", "soo_list", "ORBITAL ORBITAL (RI)") 859 ELSEIF (rigpw) THEN 860 ALLOCATE (ri_present(nkind), ri_radius(nkind)) 861 ri_present = .FALSE. 862 ri_radius = 0.0_dp 863 DO ikind = 1, nkind 864 CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis_set, basis_type="RI_HXC") 865 IF (ASSOCIATED(ri_basis_set)) THEN 866 ri_present(ikind) = .TRUE. 867 CALL get_gto_basis_set(gto_basis_set=ri_basis_set, kind_radius=ri_radius(ikind)) 868 ELSE 869 ri_present(ikind) = .FALSE. 870 ENDIF 871 END DO 872 ! set neighborlists in lri_env environment 873 CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius) 874 soo_list => qs_env%lri_env%soo_list 875 CALL build_neighbor_lists(soo_list, particle_set, atom2d, cell, pair_radius, & 876 mic=mic, molecular=molecule_only, subcells=subcells, nlname="soo_list") 877 qs_env%lri_env%soo_list => soo_list 878 ! 879 CALL pair_radius_setup(ri_present, ri_present, ri_radius, ri_radius, pair_radius) 880 saa_list => qs_env%lri_env%saa_list 881 CALL build_neighbor_lists(saa_list, particle_set, atom2d, cell, pair_radius, & 882 mic=mic, molecular=molecule_only, subcells=subcells, nlname="saa_list") 883 qs_env%lri_env%saa_list => saa_list 884 ! 885 CALL pair_radius_setup(ri_present, orb_present, ri_radius, orb_radius, pair_radius) 886 soa_list => qs_env%lri_env%soa_list 887 CALL build_neighbor_lists(soa_list, particle_set, atom2d, cell, pair_radius, & 888 mic=mic, symmetric=.FALSE., molecular=molecule_only, & 889 subcells=subcells, operator_type="ABC", nlname="saa_list") 890 qs_env%lri_env%soa_list => soa_list 891 END IF 892 893 ! Build the neighbor lists for the ALMO delocalization 894 IF (almo) THEN 895 DO ikind = 1, nkind 896 CALL get_atomic_kind(atomic_kind_set(ikind), rcov=almo_rcov, rvdw=almo_rvdw) 897 ! multiply the radius by some hard-coded number 898 c_radius(ikind) = MAX(almo_rcov, almo_rvdw)*bohr* & 899 almo_max_cutoff_multiplier 900 END DO 901 default_present = .TRUE. !include all atoms (even without basis) 902 CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius) 903 CALL build_neighbor_lists(sab_almo, particle_set, atom2d, cell, pair_radius, & 904 subcells=subcells, operator_type="PP", nlname="sab_almo") 905 CALL set_ks_env(ks_env=ks_env, sab_almo=sab_almo) 906 ENDIF 907 908 ! Print particle distribution 909 print_key_path = "PRINT%DISTRIBUTION" 910 IF (BTEST(cp_print_key_should_output(logger%iter_info, force_env_section, & 911 print_key_path), & 912 cp_p_file)) THEN 913 iw = cp_print_key_unit_nr(logger=logger, & 914 basis_section=force_env_section, & 915 print_key_path=print_key_path, & 916 extension=".out") 917 CALL write_neighbor_distribution(sab_orb, qs_kind_set, iw, para_env) 918 CALL cp_print_key_finished_output(unit_nr=iw, & 919 logger=logger, & 920 basis_section=force_env_section, & 921 print_key_path=print_key_path) 922 END IF 923 924 ! Release work storage 925 CALL atom2d_cleanup(atom2d) 926 927 DEALLOCATE (atom2d) 928 DEALLOCATE (orb_present, default_present, core_present) 929 DEALLOCATE (orb_radius, aux_fit_radius, c_radius, core_radius) 930 DEALLOCATE (calpha, zeff) 931 DEALLOCATE (pair_radius) 932 IF (gth_potential_present .OR. sgp_potential_present) THEN 933 DEALLOCATE (ppl_present, ppl_radius) 934 DEALLOCATE (ppnl_present, ppnl_radius) 935 END IF 936 IF (paw_atom_present) THEN 937 DEALLOCATE (oce_present, oce_radius) 938 ENDIF 939 IF (all_potential_present .OR. sgp_potential_present) THEN 940 DEALLOCATE (all_present, all_pot_rad) 941 END IF 942 943 CALL timestop(handle) 944 945 END SUBROUTINE build_qs_neighbor_lists 946 947! ************************************************************************************************** 948!> \brief Build simple pair neighbor lists. 949!> \param ab_list ... 950!> \param particle_set ... 951!> \param atom ... 952!> \param cell ... 953!> \param pair_radius ... 954!> \param subcells ... 955!> \param mic ... 956!> \param symmetric ... 957!> \param molecular ... 958!> \param subset_of_mol ... 959!> \param current_subset ... 960!> \param operator_type ... 961!> \param nlname ... 962!> \param atomb_to_keep the list of atom indices to keep for pairs from the atom2d%b_list 963!> \date 20.03.2002 964!> \par History 965!> - Major refactoring (25.07.2010,jhu) 966!> - Added option to filter out atoms from list_b (08.2018, A. Bussy) 967!> \author MK 968!> \version 2.0 969! ************************************************************************************************** 970 SUBROUTINE build_neighbor_lists(ab_list, particle_set, atom, cell, pair_radius, subcells, & 971 mic, symmetric, molecular, subset_of_mol, current_subset, & 972 operator_type, nlname, atomb_to_keep) 973 974 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 975 POINTER :: ab_list 976 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 977 TYPE(local_atoms_type), DIMENSION(:), INTENT(IN) :: atom 978 TYPE(cell_type), POINTER :: cell 979 REAL(dp), DIMENSION(:, :), INTENT(IN) :: pair_radius 980 REAL(dp), INTENT(IN) :: subcells 981 LOGICAL, INTENT(IN), OPTIONAL :: mic, symmetric, molecular 982 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: subset_of_mol 983 INTEGER, OPTIONAL :: current_subset 984 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: operator_type 985 CHARACTER(LEN=*), INTENT(IN) :: nlname 986 INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: atomb_to_keep 987 988 CHARACTER(len=*), PARAMETER :: routineN = 'build_neighbor_lists', & 989 routineP = moduleN//':'//routineN 990 991 TYPE local_lists 992 INTEGER, DIMENSION(:), POINTER :: list 993 END TYPE local_lists 994 INTEGER :: atom_a, atom_b, handle, i, iab, iatom, iatom_local, & 995 iatom_subcell, icell, ikind, j, jatom, jatom_local, jcell, jkind, k, & 996 kcell, maxat, mol_a, mol_b, nkind, otype, natom 997 INTEGER, DIMENSION(3) :: cell_b, ncell, & 998 nsubcell, periodic 999 INTEGER, DIMENSION(:), POINTER :: index_list 1000 LOGICAL :: include_ab, my_mic, & 1001 my_molecular, my_symmetric, my_sort_atomb 1002 LOGICAL, ALLOCATABLE, DIMENSION(:) :: pres_a, pres_b 1003 REAL(dp) :: rab2, rab2_max, rab_max, rabm, deth, subcell_scale 1004 REAL(dp), DIMENSION(3) :: r, rab, ra, rb, sab_max, sb, & 1005 sb_pbc, sb_min, sb_max, rab_pbc, pd, sab_max_guard 1006 INTEGER, ALLOCATABLE, DIMENSION(:) :: nlista, nlistb 1007 TYPE(local_lists), DIMENSION(:), POINTER :: lista, listb 1008 TYPE(neighbor_list_p_type), & 1009 ALLOCATABLE, DIMENSION(:) :: kind_a 1010 TYPE(neighbor_list_set_type), POINTER :: neighbor_list_set 1011 TYPE(subcell_type), DIMENSION(:, :, :), & 1012 POINTER :: subcell 1013 REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: r_pbc 1014 1015 CALL timeset(routineN//"_"//TRIM(nlname), handle) 1016 1017 ! input options 1018 my_mic = .FALSE. 1019 IF (PRESENT(mic)) my_mic = mic 1020 my_symmetric = .TRUE. 1021 IF (PRESENT(symmetric)) my_symmetric = symmetric 1022 my_molecular = .FALSE. 1023 ! if we have a molecular NL, MIC has to be used 1024 IF (PRESENT(molecular)) my_molecular = molecular 1025 ! check for operator types 1026 IF (PRESENT(operator_type)) THEN 1027 SELECT CASE (operator_type) 1028 CASE ("AB") 1029 otype = 1 ! simple overlap 1030 CASE ("ABC") 1031 otype = 2 ! for three center operators 1032 CPASSERT(.NOT. my_molecular) 1033 my_symmetric = .FALSE. 1034 CASE ("ABBA") 1035 otype = 3 ! for separable nonlocal operators 1036 my_symmetric = .FALSE. 1037 CASE ("PP") 1038 otype = 4 ! simple atomic pair potential list 1039 CASE default 1040 CPABORT("") 1041 END SELECT 1042 ELSE 1043 ! default is a simple AB neighbor list 1044 otype = 1 1045 END IF 1046 my_sort_atomb = .FALSE. 1047 IF (PRESENT(atomb_to_keep)) THEN 1048 my_sort_atomb = .TRUE. 1049 END IF 1050 1051 ! Deallocate the old neighbor list structure 1052 IF (ASSOCIATED(ab_list)) THEN 1053 DO iab = 1, SIZE(ab_list) 1054 CALL deallocate_neighbor_list_set(ab_list(iab)%neighbor_list_set) 1055 END DO 1056 DEALLOCATE (ab_list) 1057 END IF 1058 nkind = SIZE(atom) 1059 ! Allocate and initialize the new neighbor list structure 1060 ALLOCATE (ab_list(nkind*nkind)) 1061 DO iab = 1, SIZE(ab_list) 1062 NULLIFY (ab_list(iab)%neighbor_list_set) 1063 END DO 1064 1065 ! Allocate and initialize the kind availability 1066 ALLOCATE (pres_a(nkind), pres_b(nkind)) 1067 DO ikind = 1, nkind 1068 pres_a(ikind) = ANY(pair_radius(ikind, :) > 0._dp) 1069 pres_b(ikind) = ANY(pair_radius(:, ikind) > 0._dp) 1070 END DO 1071 1072 ! create a copy of the pbc'ed coordinates 1073 natom = SIZE(particle_set) 1074 ALLOCATE (r_pbc(3, natom)) 1075 DO i = 1, natom 1076 r_pbc(1:3, i) = pbc(particle_set(i)%r(1:3), cell) 1077 ENDDO 1078 1079 ! setup the local lists of atoms 1080 maxat = 0 1081 DO ikind = 1, nkind 1082 maxat = MAX(maxat, SIZE(atom(ikind)%list)) 1083 END DO 1084 ALLOCATE (index_list(maxat)) 1085 DO i = 1, maxat 1086 index_list(i) = i 1087 END DO 1088 ALLOCATE (lista(nkind), listb(nkind), nlista(nkind), nlistb(nkind)) 1089 nlista = 0 1090 nlistb = 0 1091 DO ikind = 1, nkind 1092 NULLIFY (lista(ikind)%list, listb(ikind)%list) 1093 SELECT CASE (otype) 1094 CASE (1) 1095 IF (ASSOCIATED(atom(ikind)%list_local_a_index)) THEN 1096 lista(ikind)%list => atom(ikind)%list_local_a_index 1097 nlista(ikind) = SIZE(lista(ikind)%list) 1098 END IF 1099 IF (ASSOCIATED(atom(ikind)%list_local_b_index)) THEN 1100 listb(ikind)%list => atom(ikind)%list_local_b_index 1101 nlistb(ikind) = SIZE(listb(ikind)%list) 1102 END IF 1103 CASE (2) 1104 IF (ASSOCIATED(atom(ikind)%list_local_a_index)) THEN 1105 lista(ikind)%list => atom(ikind)%list_local_a_index 1106 nlista(ikind) = SIZE(lista(ikind)%list) 1107 END IF 1108 nlistb(ikind) = SIZE(atom(ikind)%list) 1109 listb(ikind)%list => index_list 1110 CASE (3) 1111 CALL combine_lists(lista(ikind)%list, nlista(ikind), ikind, atom) 1112 nlistb(ikind) = SIZE(atom(ikind)%list) 1113 listb(ikind)%list => index_list 1114 CASE (4) 1115 nlista(ikind) = SIZE(atom(ikind)%list_1d) 1116 lista(ikind)%list => atom(ikind)%list_1d 1117 nlistb(ikind) = SIZE(atom(ikind)%list) 1118 listb(ikind)%list => index_list 1119 CASE default 1120 CPABORT("") 1121 END SELECT 1122 END DO 1123 1124 ! Determine max. number of local atoms 1125 maxat = 0 1126 DO ikind = 1, nkind 1127 maxat = MAX(maxat, nlista(ikind), nlistb(ikind)) 1128 END DO 1129 ALLOCATE (kind_a(2*maxat)) 1130 1131 ! Load informations about the simulation cell 1132 CALL get_cell(cell=cell, periodic=periodic, deth=deth) 1133 1134 ! Loop over all atomic kind pairs 1135 DO ikind = 1, nkind 1136 IF (.NOT. pres_a(ikind)) CYCLE 1137 1138 DO jkind = 1, nkind 1139 IF (.NOT. pres_b(jkind)) CYCLE 1140 1141 iab = ikind + nkind*(jkind - 1) 1142 1143 ! Calculate the square of the maximum interaction distance 1144 IF (pair_radius(ikind, jkind) <= 0._dp) CYCLE 1145 rab_max = pair_radius(ikind, jkind) 1146 IF (otype == 3) THEN 1147 ! Calculate the square of the maximum interaction distance 1148 ! for sac_max / ncell this must be the maximum over all kinds 1149 ! to be correct for three center terms involving different kinds 1150 rabm = MAXVAL(pair_radius(:, jkind)) 1151 ELSE 1152 rabm = rab_max 1153 END IF 1154 rab2_max = rabm*rabm 1155 1156 pd(1) = plane_distance(1, 0, 0, cell) 1157 pd(2) = plane_distance(0, 1, 0, cell) 1158 pd(3) = plane_distance(0, 0, 1, cell) 1159 1160 sab_max = rabm/pd 1161 sab_max_guard = 15.0_dp/pd 1162 1163 ! It makes sense to have fewer subcells for larger systems 1164 subcell_scale = ((125.0_dp**3)/deth)**(1.0_dp/6.0_dp) 1165 1166 ! guess the number of subcells for optimal performance, 1167 ! guard against crazy stuff triggered by very small rabm 1168 nsubcell(:) = INT(MAX(1.0_dp, MIN(0.5_dp*subcells*subcell_scale/sab_max(:), & 1169 0.5_dp*subcells*subcell_scale/sab_max_guard(:)))) 1170 1171 ! number of image cells to be considered 1172 ncell(:) = (INT(sab_max(:)) + 1)*periodic(:) 1173 1174 CALL allocate_neighbor_list_set(neighbor_list_set=ab_list(iab)%neighbor_list_set, & 1175 symmetric=my_symmetric) 1176 neighbor_list_set => ab_list(iab)%neighbor_list_set 1177 1178 DO iatom_local = 1, nlista(ikind) 1179 iatom = lista(ikind)%list(iatom_local) 1180 atom_a = atom(ikind)%list(iatom) 1181 CALL add_neighbor_list(neighbor_list_set=neighbor_list_set, & 1182 atom=atom_a, & 1183 neighbor_list=kind_a(iatom_local)%neighbor_list) 1184 END DO 1185 1186 CALL allocate_subcell(subcell, nsubcell) 1187 DO iatom_local = 1, nlista(ikind) 1188 iatom = lista(ikind)%list(iatom_local) 1189 atom_a = atom(ikind)%list(iatom) 1190 r = r_pbc(:, atom_a) 1191 CALL give_ijk_subcell(r, i, j, k, cell, nsubcell) 1192 subcell(i, j, k)%natom = subcell(i, j, k)%natom + 1 1193 END DO 1194 DO k = 1, nsubcell(3) 1195 DO j = 1, nsubcell(2) 1196 DO i = 1, nsubcell(1) 1197 maxat = subcell(i, j, k)%natom + subcell(i, j, k)%natom/10 1198 ALLOCATE (subcell(i, j, k)%atom_list(maxat)) 1199 subcell(i, j, k)%natom = 0 1200 END DO 1201 END DO 1202 END DO 1203 DO iatom_local = 1, nlista(ikind) 1204 iatom = lista(ikind)%list(iatom_local) 1205 atom_a = atom(ikind)%list(iatom) 1206 r = r_pbc(:, atom_a) 1207 CALL give_ijk_subcell(r, i, j, k, cell, nsubcell) 1208 subcell(i, j, k)%natom = subcell(i, j, k)%natom + 1 1209 subcell(i, j, k)%atom_list(subcell(i, j, k)%natom) = iatom_local 1210 END DO 1211 1212 DO jatom_local = 1, nlistb(jkind) 1213 jatom = listb(jkind)%list(jatom_local) 1214 atom_b = atom(jkind)%list(jatom) 1215 IF (my_sort_atomb .AND. .NOT. my_symmetric) THEN 1216 IF (.NOT. ANY(atomb_to_keep == atom_b)) CYCLE 1217 END IF 1218 IF (my_molecular) THEN 1219 mol_b = atom(jkind)%list_b_mol(jatom_local) 1220 IF (PRESENT(subset_of_mol)) THEN 1221 IF (subset_of_mol(mol_b) .NE. current_subset) CYCLE 1222 END IF 1223 END IF 1224 r = r_pbc(:, atom_b) 1225 CALL real_to_scaled(sb_pbc(:), r(:), cell) 1226 1227 loop2_kcell: DO kcell = -ncell(3), ncell(3) 1228 sb(3) = sb_pbc(3) + REAL(kcell, dp) 1229 sb_min(3) = sb(3) - sab_max(3) 1230 sb_max(3) = sb(3) + sab_max(3) 1231 IF (periodic(3) /= 0) THEN 1232 IF (sb_min(3) >= 0.5_dp) EXIT loop2_kcell 1233 IF (sb_max(3) < -0.5_dp) CYCLE loop2_kcell 1234 END IF 1235 cell_b(3) = kcell 1236 1237 loop2_jcell: DO jcell = -ncell(2), ncell(2) 1238 sb(2) = sb_pbc(2) + REAL(jcell, dp) 1239 sb_min(2) = sb(2) - sab_max(2) 1240 sb_max(2) = sb(2) + sab_max(2) 1241 IF (periodic(2) /= 0) THEN 1242 IF (sb_min(2) >= 0.5_dp) EXIT loop2_jcell 1243 IF (sb_max(2) < -0.5_dp) CYCLE loop2_jcell 1244 END IF 1245 cell_b(2) = jcell 1246 1247 loop2_icell: DO icell = -ncell(1), ncell(1) 1248 sb(1) = sb_pbc(1) + REAL(icell, dp) 1249 sb_min(1) = sb(1) - sab_max(1) 1250 sb_max(1) = sb(1) + sab_max(1) 1251 IF (periodic(1) /= 0) THEN 1252 IF (sb_min(1) >= 0.5_dp) EXIT loop2_icell 1253 IF (sb_max(1) < -0.5_dp) CYCLE loop2_icell 1254 END IF 1255 cell_b(1) = icell 1256 1257 CALL scaled_to_real(rb, sb, cell) 1258 1259 loop_k: DO k = 1, nsubcell(3) 1260 loop_j: DO j = 1, nsubcell(2) 1261 loop_i: DO i = 1, nsubcell(1) 1262 1263 ! FIXME for non-periodic systems, the whole subcell trick is skipped 1264 ! yielding a Natom**2 pair list build. 1265 IF (periodic(3) /= 0) THEN 1266 IF (sb_max(3) < subcell(i, j, k)%s_min(3)) EXIT loop_k 1267 IF (sb_min(3) >= subcell(i, j, k)%s_max(3)) CYCLE loop_k 1268 END IF 1269 1270 IF (periodic(2) /= 0) THEN 1271 IF (sb_max(2) < subcell(i, j, k)%s_min(2)) EXIT loop_j 1272 IF (sb_min(2) >= subcell(i, j, k)%s_max(2)) CYCLE loop_j 1273 END IF 1274 1275 IF (periodic(1) /= 0) THEN 1276 IF (sb_max(1) < subcell(i, j, k)%s_min(1)) EXIT loop_i 1277 IF (sb_min(1) >= subcell(i, j, k)%s_max(1)) CYCLE loop_i 1278 END IF 1279 1280 IF (subcell(i, j, k)%natom == 0) CYCLE 1281 1282 DO iatom_subcell = 1, subcell(i, j, k)%natom 1283 iatom_local = subcell(i, j, k)%atom_list(iatom_subcell) 1284 iatom = lista(ikind)%list(iatom_local) 1285 atom_a = atom(ikind)%list(iatom) 1286 IF (my_molecular) THEN 1287 mol_a = atom(ikind)%list_a_mol(iatom_local) 1288 IF (mol_a /= mol_b) CYCLE 1289 END IF 1290 IF (my_symmetric) THEN 1291 IF (atom_a > atom_b) THEN 1292 include_ab = (MODULO(atom_a + atom_b, 2) /= 0) 1293 ELSE 1294 include_ab = (MODULO(atom_a + atom_b, 2) == 0) 1295 END IF 1296 IF (my_sort_atomb) THEN 1297 IF ((.NOT. ANY(atomb_to_keep == atom_b)) .AND. & 1298 (.NOT. ANY(atomb_to_keep == atom_a))) THEN 1299 include_ab = .FALSE. 1300 END IF 1301 END IF 1302 ELSE 1303 include_ab = .TRUE. 1304 END IF 1305 IF (include_ab) THEN 1306 ra(:) = r_pbc(:, atom_a) 1307 rab(:) = rb(:) - ra(:) 1308 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 1309 IF (rab2 < rab2_max) THEN 1310 include_ab = .TRUE. 1311 IF (my_mic) THEN 1312 ! only if rab is minimum image the pair will be included 1313 ! ideally the range of the pair list is < L/2 so that this never triggers 1314 rab_pbc(:) = pbc(rab(:), cell) 1315 IF (SUM((rab_pbc - rab)**2) > EPSILON(1.0_dp)) THEN 1316 include_ab = .FALSE. 1317 ENDIF 1318 ENDIF 1319 IF (include_ab) THEN 1320 CALL add_neighbor_node( & 1321 neighbor_list=kind_a(iatom_local)%neighbor_list, & 1322 neighbor=atom_b, & 1323 cell=cell_b, & 1324 r=rab, & 1325 nkind=nkind) 1326 ENDIF 1327 END IF 1328 END IF 1329 END DO 1330 1331 END DO loop_i 1332 END DO loop_j 1333 END DO loop_k 1334 1335 END DO loop2_icell 1336 END DO loop2_jcell 1337 END DO loop2_kcell 1338 1339 END DO 1340 1341 CALL deallocate_subcell(subcell) 1342 1343 END DO 1344 END DO 1345 1346 SELECT CASE (otype) 1347 CASE (1:2, 4) 1348 CASE (3) 1349 DO ikind = 1, nkind 1350 DEALLOCATE (lista(ikind)%list) 1351 END DO 1352 CASE default 1353 CPABORT("") 1354 END SELECT 1355 DEALLOCATE (kind_a, pres_a, pres_b, lista, listb, nlista, nlistb) 1356 DEALLOCATE (index_list) 1357 DEALLOCATE (r_pbc) 1358 1359 CALL timestop(handle) 1360 1361 END SUBROUTINE build_neighbor_lists 1362 1363! ************************************************************************************************** 1364!> \brief Build a neighborlist 1365!> \param ab_list ... 1366!> \param basis_set_a ... 1367!> \param basis_set_b ... 1368!> \param qs_env ... 1369!> \param mic ... 1370!> \param symmetric ... 1371!> \param molecular ... 1372!> \param operator_type ... 1373!> \date 14.03.2016 1374!> \author JGH 1375! ************************************************************************************************** 1376 SUBROUTINE setup_neighbor_list(ab_list, basis_set_a, basis_set_b, qs_env, & 1377 mic, symmetric, molecular, operator_type) 1378 1379 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 1380 POINTER :: ab_list 1381 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_a 1382 TYPE(gto_basis_set_p_type), DIMENSION(:), & 1383 OPTIONAL, POINTER :: basis_set_b 1384 TYPE(qs_environment_type), POINTER :: qs_env 1385 LOGICAL, INTENT(IN), OPTIONAL :: mic, symmetric, molecular 1386 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: operator_type 1387 1388 CHARACTER(len=*), PARAMETER :: routineN = 'setup_neighbor_list', & 1389 routineP = moduleN//':'//routineN 1390 1391 CHARACTER(LEN=4) :: otype 1392 INTEGER :: ikind, nkind 1393 LOGICAL :: my_mic, my_molecular, my_symmetric 1394 LOGICAL, ALLOCATABLE, DIMENSION(:) :: a_present, b_present 1395 REAL(dp), ALLOCATABLE, DIMENSION(:) :: a_radius, b_radius 1396 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius 1397 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1398 TYPE(cell_type), POINTER :: cell 1399 TYPE(distribution_1d_type), POINTER :: distribution_1d 1400 TYPE(distribution_2d_type), POINTER :: distribution_2d 1401 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_a, basis_b 1402 TYPE(gto_basis_set_type), POINTER :: abas, bbas 1403 TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d 1404 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 1405 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1406 1407 basis_a => basis_set_a 1408 IF (PRESENT(basis_set_b)) THEN 1409 basis_b => basis_set_b 1410 my_symmetric = .FALSE. 1411 ELSE 1412 basis_b => basis_set_a 1413 my_symmetric = .TRUE. 1414 END IF 1415 IF (PRESENT(symmetric)) my_symmetric = symmetric 1416 1417 IF (PRESENT(mic)) THEN 1418 my_mic = mic 1419 ELSE 1420 my_mic = .FALSE. 1421 END IF 1422 1423 IF (PRESENT(molecular)) THEN 1424 my_molecular = molecular 1425 ELSE 1426 my_molecular = .FALSE. 1427 END IF 1428 1429 IF (PRESENT(operator_type)) THEN 1430 otype = operator_type 1431 ELSE 1432 ! default is a simple AB neighbor list 1433 otype = "AB" 1434 END IF 1435 1436 nkind = SIZE(basis_a) 1437 ALLOCATE (a_present(nkind), b_present(nkind)) 1438 a_present = .FALSE. 1439 b_present = .FALSE. 1440 ALLOCATE (a_radius(nkind), b_radius(nkind)) 1441 a_radius = 0.0_dp 1442 b_radius = 0.0_dp 1443 DO ikind = 1, nkind 1444 IF (ASSOCIATED(basis_a(ikind)%gto_basis_set)) THEN 1445 a_present(ikind) = .TRUE. 1446 abas => basis_a(ikind)%gto_basis_set 1447 CALL get_gto_basis_set(gto_basis_set=abas, kind_radius=a_radius(ikind)) 1448 END IF 1449 IF (ASSOCIATED(basis_b(ikind)%gto_basis_set)) THEN 1450 b_present(ikind) = .TRUE. 1451 bbas => basis_b(ikind)%gto_basis_set 1452 CALL get_gto_basis_set(gto_basis_set=bbas, kind_radius=b_radius(ikind)) 1453 END IF 1454 END DO 1455 1456 ALLOCATE (pair_radius(nkind, nkind)) 1457 pair_radius = 0.0_dp 1458 CALL pair_radius_setup(a_present, b_present, a_radius, b_radius, pair_radius) 1459 1460 CALL get_qs_env(qs_env, & 1461 atomic_kind_set=atomic_kind_set, & 1462 cell=cell, & 1463 distribution_2d=distribution_2d, & 1464 local_particles=distribution_1d, & 1465 particle_set=particle_set, & 1466 molecule_set=molecule_set) 1467 1468 ALLOCATE (atom2d(nkind)) 1469 CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, & 1470 molecule_set, my_molecular, particle_set=particle_set) 1471 CALL build_neighbor_lists(ab_list, particle_set, atom2d, cell, pair_radius, & 1472 mic=my_mic, symmetric=my_symmetric, molecular=my_molecular, & 1473 subcells=2.0_dp, nlname="AUX_NL") 1474 1475 CALL atom2d_cleanup(atom2d) 1476 1477 DEALLOCATE (a_present, b_present, a_radius, b_radius, pair_radius, atom2d) 1478 1479 END SUBROUTINE setup_neighbor_list 1480 1481! ************************************************************************************************** 1482!> \brief ... 1483!> \param list ... 1484!> \param n ... 1485!> \param ikind ... 1486!> \param atom ... 1487! ************************************************************************************************** 1488 SUBROUTINE combine_lists(list, n, ikind, atom) 1489 INTEGER, DIMENSION(:), POINTER :: list 1490 INTEGER, INTENT(OUT) :: n 1491 INTEGER, INTENT(IN) :: ikind 1492 TYPE(local_atoms_type), DIMENSION(:), INTENT(IN) :: atom 1493 1494 CHARACTER(len=*), PARAMETER :: routineN = 'combine_lists', routineP = moduleN//':'//routineN 1495 1496 INTEGER :: i, ib, na, nb 1497 INTEGER, DIMENSION(:), POINTER :: lista, listb 1498 1499 CPASSERT(.NOT. ASSOCIATED(list)) 1500 1501 lista => atom(ikind)%list_local_a_index 1502 listb => atom(ikind)%list_local_b_index 1503 1504 IF (ASSOCIATED(lista)) THEN 1505 na = SIZE(lista) 1506 ELSE 1507 na = 0 1508 END IF 1509 1510 IF (ASSOCIATED(listb)) THEN 1511 nb = SIZE(listb) 1512 ELSE 1513 nb = 0 1514 END IF 1515 1516 ALLOCATE (list(na + nb)) 1517 1518 n = na 1519 IF (na .GT. 0) list(1:na) = lista(1:na) 1520 IF (nb .GT. 0) THEN 1521 loopb: DO ib = 1, nb 1522 DO i = 1, na 1523 IF (listb(ib) == list(i)) CYCLE loopb 1524 END DO 1525 n = n + 1 1526 list(n) = listb(ib) 1527 END DO loopb 1528 ENDIF 1529 END SUBROUTINE combine_lists 1530 1531! ************************************************************************************************** 1532 1533! ************************************************************************************************** 1534!> \brief ... 1535!> \param present_a ... 1536!> \param present_b ... 1537!> \param radius_a ... 1538!> \param radius_b ... 1539!> \param pair_radius ... 1540! ************************************************************************************************** 1541 SUBROUTINE pair_radius_setup(present_a, present_b, radius_a, radius_b, pair_radius) 1542 LOGICAL, DIMENSION(:), INTENT(IN) :: present_a, present_b 1543 REAL(dp), DIMENSION(:), INTENT(IN) :: radius_a, radius_b 1544 REAL(dp), DIMENSION(:, :), INTENT(OUT) :: pair_radius 1545 1546 CHARACTER(len=*), PARAMETER :: routineN = 'pair_radius_setup', & 1547 routineP = moduleN//':'//routineN 1548 1549 INTEGER :: i, j, nkind 1550 1551 nkind = SIZE(present_a) 1552 1553 pair_radius = 0._dp 1554 1555 DO i = 1, nkind 1556 IF (.NOT. present_a(i)) CYCLE 1557 DO j = 1, nkind 1558 IF (.NOT. present_b(j)) CYCLE 1559 pair_radius(i, j) = radius_a(i) + radius_b(j) 1560 END DO 1561 END DO 1562 1563 END SUBROUTINE pair_radius_setup 1564 1565! ************************************************************************************************** 1566!> \brief Print the distribution of the simple pair neighbor list. 1567!> \param ab ... 1568!> \param qs_kind_set ... 1569!> \param output_unit ... 1570!> \param para_env ... 1571!> \date 19.06.2003 1572!> \author MK 1573!> \version 1.0 1574! ************************************************************************************************** 1575 SUBROUTINE write_neighbor_distribution(ab, qs_kind_set, output_unit, para_env) 1576 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 1577 POINTER :: ab 1578 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1579 INTEGER, INTENT(in) :: output_unit 1580 TYPE(cp_para_env_type), POINTER :: para_env 1581 1582 CHARACTER(len=*), PARAMETER :: routineN = 'write_neighbor_distribution', & 1583 routineP = moduleN//':'//routineN 1584 LOGICAL, PARAMETER :: full_output = .FALSE. 1585 1586 INTEGER :: group, handle, ikind, inode, ipe, jkind, & 1587 mype, n, nkind, nnode, npe 1588 INTEGER(int_8) :: nblock_max, nblock_sum, nelement_max, & 1589 nelement_sum, tmp(2) 1590 INTEGER, ALLOCATABLE, DIMENSION(:) :: nblock, nelement, nnsgf 1591 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 1592 TYPE(neighbor_list_iterator_p_type), & 1593 DIMENSION(:), POINTER :: nl_iterator 1594 1595 CALL timeset(routineN, handle) 1596 group = para_env%group 1597 mype = para_env%mepos + 1 1598 npe = para_env%num_pe 1599 1600 ! Allocate work storage 1601 ALLOCATE (nblock(npe), nelement(npe)) 1602 nblock(:) = 0 1603 nelement(:) = 0 1604 nkind = SIZE(qs_kind_set) 1605 ALLOCATE (nnsgf(nkind)) 1606 nnsgf = 1 1607 DO ikind = 1, nkind 1608 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) 1609 IF (ASSOCIATED(orb_basis_set)) THEN 1610 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, nsgf=nnsgf(ikind)) 1611 END IF 1612 END DO 1613 1614 CALL neighbor_list_iterator_create(nl_iterator, ab) 1615 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 1616 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, nnode=nnode) 1617 IF (inode == 1) THEN 1618 n = nnsgf(ikind)*nnsgf(jkind) 1619 nblock(mype) = nblock(mype) + nnode 1620 nelement(mype) = nelement(mype) + n*nnode 1621 END IF 1622 END DO 1623 CALL neighbor_list_iterator_release(nl_iterator) 1624 1625 IF (full_output) THEN 1626 ! XXXXXXXX should gather/scatter this on ionode 1627 CALL mp_sum(nblock, group) 1628 CALL mp_sum(nelement, group) 1629 1630 nblock_sum = SUM(INT(nblock, KIND=int_8)) 1631 nelement_sum = SUM(INT(nelement, KIND=int_8)) 1632 ELSE 1633 nblock_sum = nblock(mype) 1634 nblock_max = nblock(mype) 1635 nelement_sum = nelement(mype) 1636 nelement_max = nelement(mype) 1637 tmp = (/nblock_sum, nelement_sum/) 1638 CALL mp_sum(tmp, group) 1639 nblock_sum = tmp(1); nelement_sum = tmp(2) 1640 tmp = (/nblock_max, nelement_max/) 1641 CALL mp_max(tmp, group) 1642 nblock_max = tmp(1); nelement_max = tmp(2) 1643 ENDIF 1644 1645 IF (output_unit > 0) THEN 1646 IF (full_output) THEN 1647 WRITE (UNIT=output_unit, & 1648 FMT="(/,/,T2,A,/,/,T3,A,/,/,(T4,I6,T27,I10,T55,I10))") & 1649 "DISTRIBUTION OF THE NEIGHBOR LISTS", & 1650 "Process Number of particle pairs Number of matrix elements", & 1651 (ipe - 1, nblock(ipe), nelement(ipe), ipe=1, npe) 1652 WRITE (UNIT=output_unit, FMT="(/,T7,A3,T27,I10,T55,I10)") & 1653 "Sum", SUM(nblock), SUM(nelement) 1654 ELSE 1655 WRITE (UNIT=output_unit, FMT="(/,T2,A)") "DISTRIBUTION OF THE NEIGHBOR LISTS" 1656 WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Total number of particle pairs:", nblock_sum 1657 WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Total number of matrix elements:", nelement_sum 1658 WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Average number of particle pairs:", (nblock_sum + npe - 1)/npe 1659 WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Maximum number of particle pairs:", nblock_max 1660 WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Average number of matrix element:", (nelement_sum + npe - 1)/npe 1661 WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Maximum number of matrix elements:", nelement_max 1662 ENDIF 1663 END IF 1664 1665 ! Release work storage 1666 1667 DEALLOCATE (nblock, nelement, nnsgf) 1668 1669 CALL timestop(handle) 1670 1671 END SUBROUTINE write_neighbor_distribution 1672 1673! ************************************************************************************************** 1674!> \brief Write a set of neighbor lists to the output unit. 1675!> \param ab ... 1676!> \param particle_set ... 1677!> \param cell ... 1678!> \param para_env ... 1679!> \param neighbor_list_section ... 1680!> \param nl_type ... 1681!> \param middle_name ... 1682!> \param nlname ... 1683!> \date 04.03.2002 1684!> \par History 1685!> - Adapted to the new parallelized neighbor list version 1686!> (26.06.2003,MK) 1687!> \author MK 1688!> \version 1.0 1689! ************************************************************************************************** 1690 SUBROUTINE write_neighbor_lists(ab, particle_set, cell, para_env, neighbor_list_section, & 1691 nl_type, middle_name, nlname) 1692 1693 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 1694 POINTER :: ab 1695 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1696 TYPE(cell_type), POINTER :: cell 1697 TYPE(cp_para_env_type), POINTER :: para_env 1698 TYPE(section_vals_type), POINTER :: neighbor_list_section 1699 CHARACTER(LEN=*), INTENT(IN) :: nl_type, middle_name, nlname 1700 1701 CHARACTER(LEN=default_string_length) :: string, unit_str 1702 INTEGER :: iatom, inode, iw, jatom, mype, & 1703 nneighbor, nnode 1704 INTEGER, DIMENSION(3) :: cell_b 1705 REAL(dp) :: dab, unit_conv 1706 REAL(dp), DIMENSION(3) :: ra, rab, rb 1707 TYPE(cp_logger_type), POINTER :: logger 1708 TYPE(neighbor_list_iterator_p_type), & 1709 DIMENSION(:), POINTER :: nl_iterator 1710 1711 NULLIFY (logger) 1712 logger => cp_get_default_logger() 1713 IF (BTEST(cp_print_key_should_output(logger%iter_info, neighbor_list_section, & 1714 TRIM(nl_type)), & 1715 cp_p_file)) THEN 1716 iw = cp_print_key_unit_nr(logger=logger, & 1717 basis_section=neighbor_list_section, & 1718 print_key_path=TRIM(nl_type), & 1719 extension=".out", & 1720 middle_name=TRIM(middle_name), & 1721 local=.TRUE., & 1722 log_filename=.FALSE., & 1723 file_position="REWIND") 1724 mype = para_env%mepos 1725 CALL section_vals_val_get(neighbor_list_section, "UNIT", c_val=unit_str) 1726 unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) 1727 1728 ! Print headline 1729 string = "" 1730 WRITE (UNIT=string, FMT="(A,I5,A)") & 1731 TRIM(nlname)//" IN "//TRIM(unit_str)//" (PROCESS", mype, ")" 1732 CALL compress(string) 1733 IF (iw > 0) WRITE (UNIT=iw, FMT="(/,/,T2,A)") TRIM(string) 1734 1735 nneighbor = 0 1736 1737 CALL neighbor_list_iterator_create(nl_iterator, ab) 1738 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 1739 CALL get_iterator_info(nl_iterator, inode=inode, nnode=nnode, & 1740 iatom=iatom, jatom=jatom, cell=cell_b, r=rab) 1741 nneighbor = nneighbor + 1 1742 ra(:) = pbc(particle_set(iatom)%r, cell) 1743 rb(:) = ra(:) + rab(:) 1744 dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)) 1745 IF (iw > 0) THEN 1746 IF (inode == 1) THEN 1747 WRITE (UNIT=iw, FMT="(/,T2,I5,3X,I6,3X,3F12.6)") & 1748 iatom, nnode, ra(1:3)*unit_conv 1749 END IF 1750 WRITE (UNIT=iw, FMT="(T10,I6,3X,3I4,3F12.6,2X,F12.6)") & 1751 jatom, cell_b(1:3), rb(1:3)*unit_conv, dab*unit_conv 1752 END IF 1753 END DO 1754 CALL neighbor_list_iterator_release(nl_iterator) 1755 1756 string = "" 1757 WRITE (UNIT=string, FMT="(A,I12,A,I12)") & 1758 "Total number of neighbor interactions for process", mype, ":", & 1759 nneighbor 1760 CALL compress(string) 1761 IF (iw > 0) WRITE (UNIT=iw, FMT="(/,T2,A)") TRIM(string) 1762 CALL cp_print_key_finished_output(unit_nr=iw, & 1763 logger=logger, & 1764 basis_section=neighbor_list_section, & 1765 print_key_path=TRIM(nl_type), & 1766 local=.TRUE.) 1767 END IF 1768 1769 END SUBROUTINE write_neighbor_lists 1770 1771END MODULE qs_neighbor_lists 1772