1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Generate the atomic neighbor lists for FIST. 8!> \par History 9!> - build and update merged (11.02.2005,MK) 10!> - bug fix for PERIODIC NONE (24.02.06,MK) 11!> - Major rewriting (light memory neighbor lists): teo and joost (05.2006) 12!> - Completely new algorithm for the neighbor lists 13!> (faster and memory lighter) (Teo 08.2006) 14!> \author MK (19.11.2002,24.07.2003) 15!> Teodoro Laino (08.2006) - MAJOR REWRITING 16! ************************************************************************************************** 17MODULE fist_neighbor_lists 18 19 USE atomic_kind_types, ONLY: atomic_kind_type,& 20 get_atomic_kind,& 21 get_atomic_kind_set 22 USE cell_types, ONLY: cell_type,& 23 get_cell,& 24 pbc,& 25 plane_distance,& 26 scaled_to_real 27 USE cp_log_handling, ONLY: cp_get_default_logger,& 28 cp_logger_get_default_unit_nr,& 29 cp_logger_type 30 USE cp_output_handling, ONLY: cp_p_file,& 31 cp_print_key_finished_output,& 32 cp_print_key_should_output,& 33 cp_print_key_unit_nr 34 USE cp_para_types, ONLY: cp_para_env_type 35 USE cp_units, ONLY: cp_unit_from_cp2k 36 USE distribution_1d_types, ONLY: distribution_1d_type 37 USE exclusion_types, ONLY: exclusion_type 38 USE fist_neighbor_list_types, ONLY: fist_neighbor_add,& 39 fist_neighbor_deallocate,& 40 fist_neighbor_init,& 41 fist_neighbor_type,& 42 neighbor_kind_pairs_type 43 USE input_section_types, ONLY: section_vals_type,& 44 section_vals_val_get 45 USE kinds, ONLY: default_string_length,& 46 dp 47 USE mathlib, ONLY: matvec_3x3 48 USE memory_utilities, ONLY: reallocate 49 USE particle_types, ONLY: particle_type 50 USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm 51 USE string_utilities, ONLY: compress 52 USE subcell_types, ONLY: allocate_subcell,& 53 deallocate_subcell,& 54 give_ijk_subcell,& 55 reorder_atoms_subcell,& 56 subcell_type 57 USE util, ONLY: sort 58#include "./base/base_uses.f90" 59 60 IMPLICIT NONE 61 62 PRIVATE 63 64 ! Global parameters (in this module) 65 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fist_neighbor_lists' 66 67 TYPE local_atoms_type 68 INTEGER, DIMENSION(:), POINTER :: list, & 69 list_local_a_index 70 END TYPE local_atoms_type 71 72 ! Public subroutines 73 PUBLIC :: build_fist_neighbor_lists 74 75CONTAINS 76 77! ************************************************************************************************** 78!> \brief ... 79!> \param atomic_kind_set ... 80!> \param particle_set ... 81!> \param local_particles ... 82!> \param cell ... 83!> \param r_max ... 84!> \param r_minsq ... 85!> \param ei_scale14 ... 86!> \param vdw_scale14 ... 87!> \param nonbonded ... 88!> \param para_env ... 89!> \param build_from_scratch ... 90!> \param geo_check ... 91!> \param mm_section ... 92!> \param full_nl ... 93!> \param exclusions ... 94!> \par History 95!> 08.2006 created [tlaino] 96!> \author Teodoro Laino 97! ************************************************************************************************** 98 SUBROUTINE build_fist_neighbor_lists(atomic_kind_set, particle_set, & 99 local_particles, cell, r_max, r_minsq, ei_scale14, vdw_scale14, & 100 nonbonded, para_env, build_from_scratch, geo_check, mm_section, & 101 full_nl, exclusions) 102 103 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 104 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 105 TYPE(distribution_1d_type), OPTIONAL, POINTER :: local_particles 106 TYPE(cell_type), POINTER :: cell 107 REAL(dp), DIMENSION(:, :), INTENT(IN) :: r_max, r_minsq 108 REAL(KIND=DP), INTENT(IN) :: ei_scale14, vdw_scale14 109 TYPE(fist_neighbor_type), POINTER :: nonbonded 110 TYPE(cp_para_env_type), POINTER :: para_env 111 LOGICAL, INTENT(IN) :: build_from_scratch, geo_check 112 TYPE(section_vals_type), POINTER :: mm_section 113 LOGICAL, DIMENSION(:, :), OPTIONAL, POINTER :: full_nl 114 TYPE(exclusion_type), DIMENSION(:), OPTIONAL :: exclusions 115 116 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_fist_neighbor_lists', & 117 routineP = moduleN//':'//routineN 118 119 CHARACTER(LEN=default_string_length) :: kind_name, print_key_path, unit_str 120 INTEGER :: atom_a, handle, iatom_local, ikind, iw, & 121 maxatom, natom_local_a, nkind, & 122 output_unit 123 LOGICAL :: present_local_particles, & 124 print_subcell_grid 125 LOGICAL, DIMENSION(:), POINTER :: skip_kind 126 LOGICAL, DIMENSION(:, :), POINTER :: my_full_nl 127 TYPE(atomic_kind_type), POINTER :: atomic_kind 128 TYPE(cp_logger_type), POINTER :: logger 129 TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom 130 131 CALL timeset(routineN, handle) 132 NULLIFY (logger) 133 logger => cp_get_default_logger() 134 135 print_subcell_grid = .FALSE. 136 output_unit = cp_print_key_unit_nr(logger, mm_section, "PRINT%SUBCELL", & 137 extension=".Log") 138 IF (output_unit > 0) print_subcell_grid = .TRUE. 139 140 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 141 maxatom=maxatom) 142 143 present_local_particles = PRESENT(local_particles) 144 145 ! if exclusions matters local particles are present. Seems like only the exclusions 146 ! for the local particles are needed, which would imply a huge memory savings for fist 147 IF (PRESENT(exclusions)) THEN 148 CPASSERT(present_local_particles) 149 ENDIF 150 151 ! Allocate work storage 152 nkind = SIZE(atomic_kind_set) 153 ALLOCATE (atom(nkind)) 154 ALLOCATE (skip_kind(nkind)) 155 ! full_nl 156 IF (PRESENT(full_nl)) THEN 157 my_full_nl => full_nl 158 ELSE 159 ALLOCATE (my_full_nl(nkind, nkind)) 160 my_full_nl = .FALSE. 161 END IF 162 ! Initialize the local data structures 163 DO ikind = 1, nkind 164 atomic_kind => atomic_kind_set(ikind) 165 NULLIFY (atom(ikind)%list) 166 NULLIFY (atom(ikind)%list_local_a_index) 167 168 CALL get_atomic_kind(atomic_kind=atomic_kind, & 169 atom_list=atom(ikind)%list, name=kind_name) 170 skip_kind(ikind) = qmmm_ff_precond_only_qm(kind_name) 171 IF (present_local_particles) THEN 172 natom_local_a = local_particles%n_el(ikind) 173 ELSE 174 natom_local_a = SIZE(atom(ikind)%list) 175 END IF 176 IF (natom_local_a > 0) THEN 177 ALLOCATE (atom(ikind)%list_local_a_index(natom_local_a)) 178 ! Build index vector for mapping 179 DO iatom_local = 1, natom_local_a 180 IF (present_local_particles) THEN 181 atom_a = local_particles%list(ikind)%array(iatom_local) 182 ELSE 183 atom_a = atom(ikind)%list(iatom_local) 184 END IF 185 atom(ikind)%list_local_a_index(iatom_local) = atom_a 186 END DO 187 188 END IF 189 END DO 190 191 IF (build_from_scratch) THEN 192 IF (ASSOCIATED(nonbonded)) THEN 193 CALL fist_neighbor_deallocate(nonbonded) 194 END IF 195 END IF 196 197 ! Build the nonbonded neighbor lists 198 CALL build_neighbor_lists(nonbonded, particle_set, atom, cell, & 199 print_subcell_grid, output_unit, r_max, r_minsq, & 200 ei_scale14, vdw_scale14, geo_check, "NONBONDED", skip_kind, & 201 my_full_nl, exclusions) 202 203 ! Sort the list according kinds for each cell 204 CALL sort_neighbor_lists(nonbonded, nkind) 205 206 print_key_path = "PRINT%NEIGHBOR_LISTS" 207 208 IF (BTEST(cp_print_key_should_output(logger%iter_info, mm_section, print_key_path), & 209 cp_p_file)) THEN 210 iw = cp_print_key_unit_nr(logger=logger, & 211 basis_section=mm_section, & 212 print_key_path=print_key_path, & 213 extension=".out", & 214 middle_name="nonbonded_nl", & 215 local=.TRUE., & 216 log_filename=.FALSE., & 217 file_position="REWIND") 218 CALL section_vals_val_get(mm_section, TRIM(print_key_path)//"%UNIT", c_val=unit_str) 219 CALL write_neighbor_lists(nonbonded, particle_set, cell, para_env, iw, & 220 "NONBONDED NEIGHBOR LISTS", unit_str) 221 CALL cp_print_key_finished_output(unit_nr=iw, & 222 logger=logger, & 223 basis_section=mm_section, & 224 print_key_path=print_key_path, & 225 local=.TRUE.) 226 END IF 227 228 ! Release work storage 229 DO ikind = 1, nkind 230 NULLIFY (atom(ikind)%list) 231 IF (ASSOCIATED(atom(ikind)%list_local_a_index)) THEN 232 DEALLOCATE (atom(ikind)%list_local_a_index) 233 END IF 234 END DO 235 IF (PRESENT(full_nl)) THEN 236 NULLIFY (my_full_nl) 237 ELSE 238 DEALLOCATE (my_full_nl) 239 END IF 240 DEALLOCATE (atom) 241 242 DEALLOCATE (skip_kind) 243 244 CALL cp_print_key_finished_output(unit_nr=output_unit, & 245 logger=logger, & 246 basis_section=mm_section, & 247 print_key_path="PRINT%SUBCELL") 248 CALL timestop(handle) 249 250 END SUBROUTINE build_fist_neighbor_lists 251 252! ************************************************************************************************** 253!> \brief ... 254!> \param nonbonded ... 255!> \param particle_set ... 256!> \param atom ... 257!> \param cell ... 258!> \param print_subcell_grid ... 259!> \param output_unit ... 260!> \param r_max ... 261!> \param r_minsq ... 262!> \param ei_scale14 ... 263!> \param vdw_scale14 ... 264!> \param geo_check ... 265!> \param name ... 266!> \param skip_kind ... 267!> \param full_nl ... 268!> \param exclusions ... 269!> \par History 270!> 08.2006 created [tlaino] 271!> \author Teodoro Laino 272! ************************************************************************************************** 273 SUBROUTINE build_neighbor_lists(nonbonded, particle_set, atom, cell, & 274 print_subcell_grid, output_unit, r_max, r_minsq, & 275 ei_scale14, vdw_scale14, geo_check, name, skip_kind, full_nl, exclusions) 276 277 TYPE(fist_neighbor_type), POINTER :: nonbonded 278 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 279 TYPE(local_atoms_type), DIMENSION(:), INTENT(IN) :: atom 280 TYPE(cell_type), POINTER :: cell 281 LOGICAL, INTENT(IN) :: print_subcell_grid 282 INTEGER, INTENT(IN) :: output_unit 283 REAL(dp), DIMENSION(:, :), INTENT(IN) :: r_max, r_minsq 284 REAL(KIND=dp), INTENT(IN) :: ei_scale14, vdw_scale14 285 LOGICAL, INTENT(IN) :: geo_check 286 CHARACTER(LEN=*), INTENT(IN) :: name 287 LOGICAL, DIMENSION(:), POINTER :: skip_kind 288 LOGICAL, DIMENSION(:, :), POINTER :: full_nl 289 TYPE(exclusion_type), DIMENSION(:), OPTIONAL :: exclusions 290 291 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_neighbor_lists', & 292 routineP = moduleN//':'//routineN 293 294 INTEGER :: a_i, a_j, a_k, atom_a, atom_b, b_i, b_j, b_k, b_pi, b_pj, b_pk, bg_i, bg_j, bg_k, & 295 handle, i, i1, iatom_local, icell, icellmap, id_kind, ii, ii_start, ij, ij_start, ik, & 296 ik_start, ikind, imap, imax_cell, invcellmap, iw, ix, j, j1, jatom_local, jcell, jkind, & 297 jx, k, kcell, kx, natom_local_a, ncellmax, nkind, nkind00, tmpdim, xdim, ydim, zdim 298 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, work 299 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cellmap 300 INTEGER, DIMENSION(3) :: isubcell, ncell, nsubcell, periodic 301 LOGICAL :: any_full, atom_order, check_spline, & 302 is_full, subcell000 303 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :) :: sphcub 304 REAL(dp) :: rab2, rab2_max, rab2_min, rab_max 305 REAL(dp), DIMENSION(3) :: abc, cell_v, cv_b, rab, rb, sab_max 306 REAL(KIND=dp) :: ic(3), icx(3), radius, vv 307 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coord 308 TYPE(neighbor_kind_pairs_type), POINTER :: inv_neighbor_kind_pair, & 309 neighbor_kind_pair 310 TYPE(subcell_type), DIMENSION(:, :, :), POINTER :: subcell_a, subcell_b 311 312 CALL timeset(routineN, handle) 313 nkind = SIZE(atom) 314 nsubcell = 1 315 isubcell = 0 316 ncell = 0 317 any_full = ANY(full_nl) 318 CALL get_cell(cell=cell, & 319 abc=abc, & 320 periodic=periodic) 321 ! Determines the number of subcells 322 DO ikind = 1, nkind 323 DO jkind = ikind, nkind 324 ! Calculate the square of the maximum interaction distance 325 rab_max = r_max(ikind, jkind) 326 IF (skip_kind(ikind) .AND. skip_kind(jkind)) CYCLE 327 nsubcell(1) = MAX(nsubcell(1), CEILING(plane_distance(1, 0, 0, cell)/rab_max)) 328 nsubcell(2) = MAX(nsubcell(2), CEILING(plane_distance(0, 1, 0, cell)/rab_max)) 329 nsubcell(3) = MAX(nsubcell(3), CEILING(plane_distance(0, 0, 1, cell)/rab_max)) 330 END DO 331 END DO 332 ! Determines the number of periodic images and the number of interacting subcells 333 DO ikind = 1, nkind 334 DO jkind = ikind, nkind 335 IF (skip_kind(ikind) .AND. skip_kind(jkind)) CYCLE 336 ! Calculate the square of the maximum interaction distance 337 rab_max = r_max(ikind, jkind) 338 sab_max(1) = rab_max/plane_distance(1, 0, 0, cell) 339 sab_max(2) = rab_max/plane_distance(0, 1, 0, cell) 340 sab_max(3) = rab_max/plane_distance(0, 0, 1, cell) 341 ncell = MAX(ncell(:), CEILING(sab_max(:)*periodic(:))) 342 isubcell = MAX(isubcell(:), CEILING(sab_max(:)*REAL(nsubcell(:), KIND=dp))) 343 END DO 344 END DO 345 CALL fist_neighbor_init(nonbonded, ncell) 346 ! Print headline 347 IF (print_subcell_grid) THEN 348 WRITE (UNIT=output_unit, FMT="(/,/,T2,A,/)") & 349 "SUBCELL GRID INFO FOR THE "//TRIM(name)//" NEIGHBOR LISTS" 350 WRITE (UNIT=output_unit, FMT="(T4,A,10X,3I10)") " NUMBER OF SUBCELLS ::", nsubcell 351 WRITE (UNIT=output_unit, FMT="(T4,A,10X,3I10)") " NUMBER OF PERIODIC IMAGES ::", ncell 352 WRITE (UNIT=output_unit, FMT="(T4,A,10X,3I10)") " NUMBER OF INTERACTING SUBCELLS ::", isubcell 353 END IF 354 ! Allocate subcells 355 CALL allocate_subcell(subcell_a, nsubcell, cell=cell) 356 CALL allocate_subcell(subcell_b, nsubcell, cell=cell) 357 ! Let's map the sequence of the periodic images 358 ncellmax = MAXVAL(ncell) 359 ALLOCATE (cellmap(-ncellmax:ncellmax, -ncellmax:ncellmax, -ncellmax:ncellmax)) 360 cellmap = -1 361 imap = 0 362 nkind00 = nkind*(nkind + 1)/2 363 DO imax_cell = 0, ncellmax 364 DO kcell = -imax_cell, imax_cell 365 DO jcell = -imax_cell, imax_cell 366 DO icell = -imax_cell, imax_cell 367 IF (cellmap(icell, jcell, kcell) == -1) THEN 368 imap = imap + 1 369 cellmap(icell, jcell, kcell) = imap 370 CPASSERT(imap <= nonbonded%nlists) 371 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(imap) 372 373 neighbor_kind_pair%cell_vector(1) = icell 374 neighbor_kind_pair%cell_vector(2) = jcell 375 neighbor_kind_pair%cell_vector(3) = kcell 376 ENDIF 377 ENDDO 378 ENDDO 379 ENDDO 380 ENDDO 381 ! Mapping the spherical interaction between subcells 382 ALLOCATE (sphcub(-isubcell(1):isubcell(1), & 383 -isubcell(2):isubcell(2), & 384 -isubcell(3):isubcell(3))) 385 sphcub = .FALSE. 386 IF (ALL(isubcell /= 0)) THEN 387 radius = REAL(isubcell(1), KIND=dp)**2 + REAL(isubcell(2), KIND=dp)**2 + & 388 REAL(isubcell(3), KIND=dp)**2 389 loop1: DO k = -isubcell(3), isubcell(3) 390 loop2: DO j = -isubcell(2), isubcell(2) 391 loop3: DO i = -isubcell(1), isubcell(1) 392 ic = REAL((/i, j, k/), KIND=dp) 393 ! subcell cube vertex 394 DO kx = -1, 1 395 icx(3) = ic(3) + SIGN(0.5_dp, REAL(kx, KIND=dp)) 396 DO jx = -1, 1 397 icx(2) = ic(2) + SIGN(0.5_dp, REAL(jx, KIND=dp)) 398 DO ix = -1, 1 399 icx(1) = ic(1) + SIGN(0.5_dp, REAL(ix, KIND=dp)) 400 vv = icx(1)*icx(1) + icx(2)*icx(2) + icx(3)*icx(3) 401 vv = vv/radius 402 IF (vv <= 1.0_dp) THEN 403 sphcub(i, j, k) = .TRUE. 404 CYCLE loop3 405 END IF 406 END DO 407 END DO 408 END DO 409 END DO loop3 410 END DO loop2 411 END DO loop1 412 END IF 413 ! Mapping locally all atoms in the zeroth cell 414 ALLOCATE (coord(3, SIZE(particle_set))) 415 DO atom_a = 1, SIZE(particle_set) 416 coord(:, atom_a) = pbc(particle_set(atom_a)%r, cell) 417 END DO 418 ! Associate particles to subcells (local particles) 419 DO ikind = 1, nkind 420 IF (.NOT. ASSOCIATED(atom(ikind)%list_local_a_index)) CYCLE 421 natom_local_a = SIZE(atom(ikind)%list_local_a_index) 422 DO iatom_local = 1, natom_local_a 423 atom_a = atom(ikind)%list_local_a_index(iatom_local) 424 CALL give_ijk_subcell(coord(:, atom_a), i, j, k, cell, nsubcell) 425 subcell_a(i, j, k)%natom = subcell_a(i, j, k)%natom + 1 426 END DO 427 END DO 428 DO k = 1, nsubcell(3) 429 DO j = 1, nsubcell(2) 430 DO i = 1, nsubcell(1) 431 ALLOCATE (subcell_a(i, j, k)%atom_list(subcell_a(i, j, k)%natom)) 432 subcell_a(i, j, k)%natom = 0 433 END DO 434 END DO 435 END DO 436 DO ikind = 1, nkind 437 IF (.NOT. ASSOCIATED(atom(ikind)%list_local_a_index)) CYCLE 438 natom_local_a = SIZE(atom(ikind)%list_local_a_index) 439 DO iatom_local = 1, natom_local_a 440 atom_a = atom(ikind)%list_local_a_index(iatom_local) 441 CALL give_ijk_subcell(coord(:, atom_a), i, j, k, cell, nsubcell) 442 subcell_a(i, j, k)%natom = subcell_a(i, j, k)%natom + 1 443 subcell_a(i, j, k)%atom_list(subcell_a(i, j, k)%natom) = atom_a 444 END DO 445 END DO 446 ! Associate particles to subcells (distributed particles) 447 DO atom_b = 1, SIZE(particle_set) 448 CALL give_ijk_subcell(coord(:, atom_b), i, j, k, cell, nsubcell) 449 subcell_b(i, j, k)%natom = subcell_b(i, j, k)%natom + 1 450 END DO 451 DO k = 1, nsubcell(3) 452 DO j = 1, nsubcell(2) 453 DO i = 1, nsubcell(1) 454 ALLOCATE (subcell_b(i, j, k)%atom_list(subcell_b(i, j, k)%natom)) 455 subcell_b(i, j, k)%natom = 0 456 END DO 457 END DO 458 END DO 459 DO atom_b = 1, SIZE(particle_set) 460 CALL give_ijk_subcell(coord(:, atom_b), i, j, k, cell, nsubcell) 461 subcell_b(i, j, k)%natom = subcell_b(i, j, k)%natom + 1 462 subcell_b(i, j, k)%atom_list(subcell_b(i, j, k)%natom) = atom_b 463 END DO 464 ! Reorder atoms associated to subcells 465 tmpdim = MAXVAL(subcell_a(:, :, :)%natom) 466 tmpdim = MAX(tmpdim, MAXVAL(subcell_b(:, :, :)%natom)) 467 ALLOCATE (work(3*tmpdim)) 468 ALLOCATE (kind_of(SIZE(particle_set))) 469 DO i = 1, SIZE(particle_set) 470 kind_of(i) = particle_set(i)%atomic_kind%kind_number 471 END DO 472 DO k = 1, nsubcell(3) 473 DO j = 1, nsubcell(2) 474 DO i = 1, nsubcell(1) 475 CALL reorder_atoms_subcell(subcell_a(i, j, k)%atom_list, kind_of, work) 476 CALL reorder_atoms_subcell(subcell_b(i, j, k)%atom_list, kind_of, work) 477 END DO 478 END DO 479 END DO 480 DEALLOCATE (work, kind_of) 481 zdim = nsubcell(3) 482 ydim = nsubcell(2) 483 xdim = nsubcell(1) 484 is_full = .FALSE. 485 ! We can skip until ik>=0.. this prescreens the order of the subcells 486 ik_start = -isubcell(3) 487 IF (.NOT. any_full) ik_start = 0 488 ! Loop over first subcell 489 loop_a_k: DO a_k = 1, nsubcell(3) 490 loop_a_j: DO a_j = 1, nsubcell(2) 491 loop_a_i: DO a_i = 1, nsubcell(1) 492 IF (subcell_a(a_i, a_j, a_k)%natom == 0) CYCLE 493 ! Loop over second subcell 494 loop_b_k: DO ik = ik_start, isubcell(3) 495 bg_k = a_k + ik 496 b_k = MOD(bg_k, zdim) 497 b_pk = bg_k/zdim 498 IF (b_k <= 0) THEN 499 b_k = zdim + b_k 500 b_pk = b_pk - 1 501 END IF 502 IF ((periodic(3) == 0) .AND. (ABS(b_pk) > 0)) CYCLE 503 ! Setup the starting point.. this prescreens the order of the subcells 504 ij_start = -isubcell(2) 505 IF ((ik == 0) .AND. (ik_start == 0)) ij_start = 0 506 loop_b_j: DO ij = ij_start, isubcell(2) 507 bg_j = a_j + ij 508 b_j = MOD(bg_j, ydim) 509 b_pj = bg_j/ydim 510 IF (b_j <= 0) THEN 511 b_j = ydim + b_j 512 b_pj = b_pj - 1 513 END IF 514 IF ((periodic(2) == 0) .AND. (ABS(b_pj) > 0)) CYCLE 515 ! Setup the starting point.. this prescreens the order of the subcells 516 ii_start = -isubcell(1) 517 IF ((ij == 0) .AND. (ij_start == 0)) ii_start = 0 518 loop_b_i: DO ii = ii_start, isubcell(1) 519 ! Ellipsoidal screening of subcells 520 IF (.NOT. sphcub(ii, ij, ik)) CYCLE 521 bg_i = a_i + ii 522 b_i = MOD(bg_i, xdim) 523 b_pi = bg_i/xdim 524 IF (b_i <= 0) THEN 525 b_i = xdim + b_i 526 b_pi = b_pi - 1 527 END IF 528 IF ((periodic(1) == 0) .AND. (ABS(b_pi) > 0)) CYCLE 529 IF (subcell_b(b_i, b_j, b_k)%natom == 0) CYCLE 530 ! Find the proper neighbor kind pair 531 icellmap = cellmap(b_pi, b_pj, b_pk) 532 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(icellmap) 533 ! Find the replica vector 534 cell_v = 0.0_dp 535 IF ((b_pi /= 0) .OR. (b_pj /= 0) .OR. (b_pk /= 0)) THEN 536 cv_b(1) = b_pi; cv_b(2) = b_pj; cv_b(3) = b_pk 537 CALL scaled_to_real(cell_v, cv_b, cell) 538 END IF 539 subcell000 = (a_k == bg_k) .AND. (a_j == bg_j) .AND. (a_i == bg_i) 540 ! Loop over particles inside subcell_a and subcell_b 541 DO jatom_local = 1, subcell_b(b_i, b_j, b_k)%natom 542 atom_b = subcell_b(b_i, b_j, b_k)%atom_list(jatom_local) 543 jkind = particle_set(atom_b)%atomic_kind%kind_number 544 rb(1) = coord(1, atom_b) + cell_v(1) 545 rb(2) = coord(2, atom_b) + cell_v(2) 546 rb(3) = coord(3, atom_b) + cell_v(3) 547 DO iatom_local = 1, subcell_a(a_i, a_j, a_k)%natom 548 atom_a = subcell_a(a_i, a_j, a_k)%atom_list(iatom_local) 549 ikind = particle_set(atom_a)%atomic_kind%kind_number 550 ! Screen interaction to avoid double counting 551 atom_order = (atom_a <= atom_b) 552 ! Special case for kind combination requiring the full NL 553 IF (any_full) THEN 554 is_full = full_nl(ikind, jkind) 555 IF (is_full) THEN 556 atom_order = (atom_a == atom_b) 557 ELSE 558 IF (ik < 0) CYCLE 559 IF (ik == 0 .AND. ij < 0) CYCLE 560 IF (ij == 0 .AND. ii < 0) CYCLE 561 END IF 562 END IF 563 IF (subcell000 .AND. atom_order) CYCLE 564 rab(1) = rb(1) - coord(1, atom_a) 565 rab(2) = rb(2) - coord(2, atom_a) 566 rab(3) = rb(3) - coord(3, atom_a) 567 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 568 rab_max = r_max(ikind, jkind) 569 rab2_max = rab_max*rab_max 570 IF (rab2 < rab2_max) THEN 571 ! Diagonal storage 572 j1 = MIN(ikind, jkind) 573 i1 = MAX(ikind, jkind) - j1 + 1 574 j1 = nkind - j1 + 1 575 id_kind = nkind00 - (j1*(j1 + 1)/2) + i1 576 ! Store the pair 577 CALL fist_neighbor_add(neighbor_kind_pair, atom_a, atom_b, & 578 rab=rab, & 579 check_spline=check_spline, id_kind=id_kind, & 580 skip=(skip_kind(ikind) .AND. skip_kind(jkind)), & 581 cell=cell, ei_scale14=ei_scale14, & 582 vdw_scale14=vdw_scale14, exclusions=exclusions) 583 ! This is to handle properly when interaction radius is larger than cell size 584 IF ((atom_a == atom_b) .AND. (ik_start == 0)) THEN 585 invcellmap = cellmap(-b_pi, -b_pj, -b_pk) 586 inv_neighbor_kind_pair => nonbonded%neighbor_kind_pairs(invcellmap) 587 rab = rab - 2.0_dp*cell_v 588 CALL fist_neighbor_add(inv_neighbor_kind_pair, atom_a, atom_b, & 589 rab=rab, & 590 check_spline=check_spline, id_kind=id_kind, & 591 skip=(skip_kind(ikind) .AND. skip_kind(jkind)), & 592 cell=cell, ei_scale14=ei_scale14, & 593 vdw_scale14=vdw_scale14, exclusions=exclusions) 594 END IF 595 ! Check for too close hits 596 IF (check_spline) THEN 597 rab2_min = r_minsq(ikind, jkind) 598 IF (rab2 < rab2_min) THEN 599 iw = cp_logger_get_default_unit_nr() 600 WRITE (iw, '(T2,A,2I7,2(A,F15.8),A)') "WARNING| Particles: ", & 601 atom_a, atom_b, & 602 " at distance [au]:", SQRT(rab2), " less than: ", & 603 SQRT(rab2_min), & 604 "; increase EMAX_SPLINE." 605 IF (rab2 < rab2_min/(1.06_dp)**2) THEN 606 IF (geo_check) THEN 607 CPABORT("GEOMETRY wrong or EMAX_SPLINE too small!") 608 END IF 609 END IF 610 END IF 611 END IF 612 END IF 613 END DO 614 END DO 615 END DO loop_b_i 616 END DO loop_b_j 617 END DO loop_b_k 618 END DO loop_a_i 619 END DO loop_a_j 620 END DO loop_a_k 621 DEALLOCATE (coord) 622 DEALLOCATE (cellmap) 623 DEALLOCATE (sphcub) 624 CALL deallocate_subcell(subcell_a) 625 CALL deallocate_subcell(subcell_b) 626 627 CALL timestop(handle) 628 END SUBROUTINE build_neighbor_lists 629 630! ************************************************************************************************** 631!> \brief Write a set of neighbor lists to the output unit. 632!> \param nonbonded ... 633!> \param particle_set ... 634!> \param cell ... 635!> \param para_env ... 636!> \param output_unit ... 637!> \param name ... 638!> \param unit_str ... 639!> \par History 640!> 08.2006 created [tlaino] 641!> \author Teodoro Laino 642! ************************************************************************************************** 643 SUBROUTINE write_neighbor_lists(nonbonded, particle_set, cell, para_env, output_unit, & 644 name, unit_str) 645 646 TYPE(fist_neighbor_type), POINTER :: nonbonded 647 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 648 TYPE(cell_type), POINTER :: cell 649 TYPE(cp_para_env_type), POINTER :: para_env 650 INTEGER, INTENT(IN) :: output_unit 651 CHARACTER(LEN=*), INTENT(IN) :: name, unit_str 652 653 CHARACTER(LEN=default_string_length) :: string 654 INTEGER :: atom_a, atom_b, iab, ilist, mype, & 655 nneighbor 656 LOGICAL :: print_headline 657 REAL(dp) :: conv, dab 658 REAL(dp), DIMENSION(3) :: cell_v, ra, rab, rb 659 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair 660 661 mype = para_env%mepos 662 ! Print headline 663 string = "" 664 WRITE (UNIT=string, FMT="(A,I5,A)") & 665 TRIM(name)//" IN "//TRIM(unit_str)//" (PROCESS", mype, ")" 666 CALL compress(string) 667 IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") TRIM(string) 668 669 print_headline = .TRUE. 670 nneighbor = 0 671 conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) 672 DO iab = 1, SIZE(nonbonded%neighbor_kind_pairs) 673 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(iab) 674 CALL matvec_3x3(cell_v, cell%hmat, REAL(neighbor_kind_pair%cell_vector, KIND=dp)) 675 DO ilist = 1, neighbor_kind_pair%npairs 676 nneighbor = nneighbor + 1 677 IF (output_unit > 0) THEN 678 ! Print second part of the headline 679 atom_a = neighbor_kind_pair%list(1, ilist) 680 atom_b = neighbor_kind_pair%list(2, ilist) 681 IF (print_headline) THEN 682 WRITE (UNIT=output_unit, FMT="(T3,2(A6,3(5X,A,5X)),1X,A11,10X,A8,A5,A10,A9)") & 683 "Atom-A", "X", "Y", "Z", "Atom-B", "X", "Y", "Z", "Cell(i,j,k)", & 684 "Distance", "ONFO", "VDW-scale", "EI-scale" 685 print_headline = .FALSE. 686 END IF 687 688 ra(:) = pbc(particle_set(atom_a)%r, cell) 689 rb(:) = pbc(particle_set(atom_b)%r, cell) 690 rab = rb(:) - ra(:) + cell_v 691 dab = SQRT(DOT_PRODUCT(rab, rab)) 692 IF (ilist <= neighbor_kind_pair%nscale) THEN 693 WRITE (UNIT=output_unit, FMT="(T3,2(I6,3(1X,F10.6)),3(1X,I3),10X,F8.4,L4,F11.5,F9.5)") & 694 atom_a, ra(1:3)*conv, & 695 atom_b, rb(1:3)*conv, & 696 neighbor_kind_pair%cell_vector, & 697 dab*conv, & 698 neighbor_kind_pair%is_onfo(ilist), & 699 neighbor_kind_pair%vdw_scale(ilist), & 700 neighbor_kind_pair%ei_scale(ilist) 701 ELSE 702 WRITE (UNIT=output_unit, FMT="(T3,2(I6,3(1X,F10.6)),3(1X,I3),10X,F8.4)") & 703 atom_a, ra(1:3)*conv, & 704 atom_b, rb(1:3)*conv, & 705 neighbor_kind_pair%cell_vector, & 706 dab*conv 707 END IF 708 END IF 709 END DO ! ilist 710 END DO ! iab 711 712 string = "" 713 WRITE (UNIT=string, FMT="(A,I12,A,I12)") & 714 "Total number of neighbor interactions for process", mype, ":", & 715 nneighbor 716 CALL compress(string) 717 IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,T2,A)") TRIM(string) 718 719 END SUBROUTINE write_neighbor_lists 720 721! ************************************************************************************************** 722!> \brief Sort the generated neighbor list according the kind 723!> \param nonbonded ... 724!> \param nkinds ... 725!> \par History 726!> 09.2007 created [tlaino] University of Zurich - Reducing memory usage 727!> for the FIST neighbor lists 728!> \author Teodoro Laino - University of Zurich 729! ************************************************************************************************** 730 SUBROUTINE sort_neighbor_lists(nonbonded, nkinds) 731 732 TYPE(fist_neighbor_type), POINTER :: nonbonded 733 INTEGER, INTENT(IN) :: nkinds 734 735 CHARACTER(LEN=*), PARAMETER :: routineN = 'sort_neighbor_lists', & 736 routineP = moduleN//':'//routineN 737 738 INTEGER :: handle, iab, id_kind, ikind, ipair, & 739 jkind, max_alloc_size, npairs, nscale, & 740 tmp 741 INTEGER, ALLOCATABLE, DIMENSION(:) :: indj 742 INTEGER, DIMENSION(:), POINTER :: work 743 INTEGER, DIMENSION(:, :), POINTER :: list_copy 744 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair 745 746 NULLIFY (neighbor_kind_pair) 747 CALL timeset(routineN, handle) 748 ! define a lookup table to get jkind for a given id_kind 749 ALLOCATE (indj(nkinds*(nkinds + 1)/2)) 750 id_kind = 0 751 DO jkind = 1, nkinds 752 DO ikind = jkind, nkinds 753 id_kind = id_kind + 1 754 indj(id_kind) = jkind 755 END DO 756 END DO 757 ! loop over all nlists and sort the pairs within each list. 758 DO iab = 1, nonbonded%nlists 759 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(iab) 760 npairs = neighbor_kind_pair%npairs 761 nscale = neighbor_kind_pair%nscale 762 IF (npairs /= 0) THEN 763 IF (npairs > nscale) THEN 764 ! 1) Sort the atom pairs by id_kind. Pairs whose interactions are 765 ! scaled (possibly to zero for exclusion) are not touched. They 766 ! stay packed in the beginning. Sorting is skipped altogether when 767 ! all pairs have scaled interactions. 768 ALLOCATE (work(1:npairs - nscale)) 769 ALLOCATE (list_copy(2, 1:npairs - nscale)) 770 ! Copy of the pair list is required to perform the permutation below 771 ! correctly. 772 list_copy = neighbor_kind_pair%list(:, nscale + 1:npairs) 773 CALL sort(neighbor_kind_pair%id_kind(nscale + 1:npairs), npairs - nscale, work) 774 ! Reorder atoms using the same permutation that was used to sort 775 ! the array id_kind. 776 DO ipair = nscale + 1, npairs 777 tmp = work(ipair - nscale) 778 neighbor_kind_pair%list(1, ipair) = list_copy(1, tmp) 779 neighbor_kind_pair%list(2, ipair) = list_copy(2, tmp) 780 END DO 781 DEALLOCATE (work) 782 DEALLOCATE (list_copy) 783 END IF 784 ! 2) determine the intervals (groups) in the pair list that correspond 785 ! to a certain id_kind. also store the corresponding ikind and 786 ! jkind. Note that this part does not assume ikind to be sorted, 787 ! but it only makes sense when contiguous blobs of the same ikind 788 ! are present. 789 ! Allocate sufficient memory in case all pairs of atom kinds are 790 ! present, and also provide storage for the pairs with exclusion 791 ! flags, which are unsorted. 792 max_alloc_size = nkinds*(nkinds + 1)/2 + nscale 793 IF (ASSOCIATED(neighbor_kind_pair%grp_kind_start)) THEN 794 DEALLOCATE (neighbor_kind_pair%grp_kind_start) 795 END IF 796 ALLOCATE (neighbor_kind_pair%grp_kind_start(max_alloc_size)) 797 IF (ASSOCIATED(neighbor_kind_pair%grp_kind_end)) THEN 798 DEALLOCATE (neighbor_kind_pair%grp_kind_end) 799 END IF 800 ALLOCATE (neighbor_kind_pair%grp_kind_end(max_alloc_size)) 801 IF (ASSOCIATED(neighbor_kind_pair%ij_kind)) THEN 802 DEALLOCATE (neighbor_kind_pair%ij_kind) 803 END IF 804 ALLOCATE (neighbor_kind_pair%ij_kind(2, max_alloc_size)) 805 ! Start the first interval. 806 ipair = 1 807 neighbor_kind_pair%ngrp_kind = 1 808 neighbor_kind_pair%grp_kind_start(neighbor_kind_pair%ngrp_kind) = ipair 809 ! Get ikind and jkind corresponding to id_kind. 810 id_kind = neighbor_kind_pair%id_kind(ipair) 811 jkind = indj(id_kind) 812 tmp = nkinds - jkind 813 ikind = nkinds + id_kind - nkinds*(nkinds + 1)/2 + (tmp*(tmp + 1)/2) 814 neighbor_kind_pair%ij_kind(1, neighbor_kind_pair%ngrp_kind) = ikind 815 neighbor_kind_pair%ij_kind(2, neighbor_kind_pair%ngrp_kind) = jkind 816 ! Define the remaining intervals. 817 DO ipair = 2, npairs 818 IF (neighbor_kind_pair%id_kind(ipair) /= neighbor_kind_pair%id_kind(ipair - 1)) THEN 819 neighbor_kind_pair%grp_kind_end(neighbor_kind_pair%ngrp_kind) = ipair - 1 820 neighbor_kind_pair%ngrp_kind = neighbor_kind_pair%ngrp_kind + 1 821 neighbor_kind_pair%grp_kind_start(neighbor_kind_pair%ngrp_kind) = ipair 822 ! Get ikind and jkind corresponding to id_kind. 823 id_kind = neighbor_kind_pair%id_kind(ipair) 824 jkind = indj(id_kind) 825 tmp = nkinds - jkind 826 ikind = nkinds + id_kind - nkinds*(nkinds + 1)/2 + (tmp*(tmp + 1)/2) 827 neighbor_kind_pair%ij_kind(1, neighbor_kind_pair%ngrp_kind) = ikind 828 neighbor_kind_pair%ij_kind(2, neighbor_kind_pair%ngrp_kind) = jkind 829 END IF 830 END DO 831 ! Finish the last interval. 832 neighbor_kind_pair%grp_kind_end(neighbor_kind_pair%ngrp_kind) = npairs 833 ! Reduce the grp arrays to the actual size because not all pairs of 834 ! atom types have to be present in this pair list. 835 CALL reallocate(neighbor_kind_pair%grp_kind_start, 1, neighbor_kind_pair%ngrp_kind) 836 CALL reallocate(neighbor_kind_pair%grp_kind_end, 1, neighbor_kind_pair%ngrp_kind) 837 CALL reallocate(neighbor_kind_pair%ij_kind, 1, 2, 1, neighbor_kind_pair%ngrp_kind) 838 END IF 839 ! Clean the memory.. 840 DEALLOCATE (neighbor_kind_pair%id_kind) 841 END DO 842 DEALLOCATE (indj) 843 CALL timestop(handle) 844 END SUBROUTINE sort_neighbor_lists 845 846END MODULE fist_neighbor_lists 847