(git:374b731)
Loading...
Searching...
No Matches
fist_neighbor_lists.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Generate the atomic neighbor lists for FIST.
10!> \par History
11!> - build and update merged (11.02.2005,MK)
12!> - bug fix for PERIODIC NONE (24.02.06,MK)
13!> - Major rewriting (light memory neighbor lists): teo and joost (05.2006)
14!> - Completely new algorithm for the neighbor lists
15!> (faster and memory lighter) (Teo 08.2006)
16!> \author MK (19.11.2002,24.07.2003)
17!> Teodoro Laino (08.2006) - MAJOR REWRITING
18! **************************************************************************************************
20
24 USE cell_types, ONLY: cell_type,&
25 get_cell,&
26 pbc,&
32 USE cp_output_handling, ONLY: cp_p_file,&
46 USE kinds, ONLY: default_string_length,&
47 dp
52 USE string_utilities, ONLY: compress
58 USE util, ONLY: sort
59#include "./base/base_uses.f90"
60
61 IMPLICIT NONE
62
63 PRIVATE
64
65 ! Global parameters (in this module)
66 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fist_neighbor_lists'
67
68 TYPE local_atoms_type
69 INTEGER, DIMENSION(:), POINTER :: list, &
70 list_local_a_index
71 END TYPE local_atoms_type
72
73 ! Public subroutines
75
76CONTAINS
77
78! **************************************************************************************************
79!> \brief ...
80!> \param atomic_kind_set ...
81!> \param particle_set ...
82!> \param local_particles ...
83!> \param cell ...
84!> \param r_max ...
85!> \param r_minsq ...
86!> \param ei_scale14 ...
87!> \param vdw_scale14 ...
88!> \param nonbonded ...
89!> \param para_env ...
90!> \param build_from_scratch ...
91!> \param geo_check ...
92!> \param mm_section ...
93!> \param full_nl ...
94!> \param exclusions ...
95!> \par History
96!> 08.2006 created [tlaino]
97!> \author Teodoro Laino
98! **************************************************************************************************
99 SUBROUTINE build_fist_neighbor_lists(atomic_kind_set, particle_set, &
100 local_particles, cell, r_max, r_minsq, ei_scale14, vdw_scale14, &
101 nonbonded, para_env, build_from_scratch, geo_check, mm_section, &
102 full_nl, exclusions)
103
104 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
105 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
106 TYPE(distribution_1d_type), OPTIONAL, POINTER :: local_particles
107 TYPE(cell_type), POINTER :: cell
108 REAL(dp), DIMENSION(:, :), INTENT(IN) :: r_max, r_minsq
109 REAL(kind=dp), INTENT(IN) :: ei_scale14, vdw_scale14
110 TYPE(fist_neighbor_type), POINTER :: nonbonded
111 TYPE(mp_para_env_type), POINTER :: para_env
112 LOGICAL, INTENT(IN) :: build_from_scratch, geo_check
113 TYPE(section_vals_type), POINTER :: mm_section
114 LOGICAL, DIMENSION(:, :), OPTIONAL, POINTER :: full_nl
115 TYPE(exclusion_type), DIMENSION(:), OPTIONAL :: exclusions
116
117 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_fist_neighbor_lists'
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 END IF
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
293 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, &
294 handle, i, i1, iatom_local, icell, icellmap, id_kind, ii, ii_start, ij, ij_start, ik, &
295 ik_start, ikind, imap, imax_cell, invcellmap, iw, ix, j, j1, jatom_local, jcell, jkind, &
296 jx, k, kcell, kx, natom_local_a, ncellmax, nkind, nkind00, tmpdim, xdim, ydim, zdim
297 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, work
298 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cellmap
299 INTEGER, DIMENSION(3) :: isubcell, ncell, nsubcell, periodic
300 LOGICAL :: any_full, atom_order, check_spline, &
301 is_full, subcell000
302 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :) :: sphcub
303 REAL(dp) :: rab2, rab2_max, rab2_min, rab_max
304 REAL(dp), DIMENSION(3) :: abc, cell_v, cv_b, rab, rb, sab_max
305 REAL(kind=dp) :: ic(3), icx(3), radius, vv
306 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coord
307 TYPE(neighbor_kind_pairs_type), POINTER :: inv_neighbor_kind_pair, &
308 neighbor_kind_pair
309 TYPE(subcell_type), DIMENSION(:, :, :), POINTER :: subcell_a, subcell_b
310
311 CALL timeset(routinen, handle)
312 nkind = SIZE(atom)
313 nsubcell = 1
314 isubcell = 0
315 ncell = 0
316 any_full = any(full_nl)
317 CALL get_cell(cell=cell, &
318 abc=abc, &
319 periodic=periodic)
320 ! Determines the number of subcells
321 DO ikind = 1, nkind
322 DO jkind = ikind, nkind
323 ! Calculate the square of the maximum interaction distance
324 rab_max = r_max(ikind, jkind)
325 IF (skip_kind(ikind) .AND. skip_kind(jkind)) cycle
326 nsubcell(1) = max(nsubcell(1), ceiling(plane_distance(1, 0, 0, cell)/rab_max))
327 nsubcell(2) = max(nsubcell(2), ceiling(plane_distance(0, 1, 0, cell)/rab_max))
328 nsubcell(3) = max(nsubcell(3), ceiling(plane_distance(0, 0, 1, cell)/rab_max))
329 END DO
330 END DO
331 ! Determines the number of periodic images and the number of interacting subcells
332 DO ikind = 1, nkind
333 DO jkind = ikind, nkind
334 IF (skip_kind(ikind) .AND. skip_kind(jkind)) cycle
335 ! Calculate the square of the maximum interaction distance
336 rab_max = r_max(ikind, jkind)
337 sab_max(1) = rab_max/plane_distance(1, 0, 0, cell)
338 sab_max(2) = rab_max/plane_distance(0, 1, 0, cell)
339 sab_max(3) = rab_max/plane_distance(0, 0, 1, cell)
340 ncell = max(ncell(:), ceiling(sab_max(:)*periodic(:)))
341 isubcell = max(isubcell(:), ceiling(sab_max(:)*real(nsubcell(:), kind=dp)))
342 END DO
343 END DO
344 CALL fist_neighbor_init(nonbonded, ncell)
345 ! Print headline
346 IF (print_subcell_grid) THEN
347 WRITE (unit=output_unit, fmt="(/,/,T2,A,/)") &
348 "SUBCELL GRID INFO FOR THE "//trim(name)//" NEIGHBOR LISTS"
349 WRITE (unit=output_unit, fmt="(T4,A,10X,3I10)") " NUMBER OF SUBCELLS ::", nsubcell
350 WRITE (unit=output_unit, fmt="(T4,A,10X,3I10)") " NUMBER OF PERIODIC IMAGES ::", ncell
351 WRITE (unit=output_unit, fmt="(T4,A,10X,3I10)") " NUMBER OF INTERACTING SUBCELLS ::", isubcell
352 END IF
353 ! Allocate subcells
354 CALL allocate_subcell(subcell_a, nsubcell, cell=cell)
355 CALL allocate_subcell(subcell_b, nsubcell, cell=cell)
356 ! Let's map the sequence of the periodic images
357 ncellmax = maxval(ncell)
358 ALLOCATE (cellmap(-ncellmax:ncellmax, -ncellmax:ncellmax, -ncellmax:ncellmax))
359 cellmap = -1
360 imap = 0
361 nkind00 = nkind*(nkind + 1)/2
362 DO imax_cell = 0, ncellmax
363 DO kcell = -imax_cell, imax_cell
364 DO jcell = -imax_cell, imax_cell
365 DO icell = -imax_cell, imax_cell
366 IF (cellmap(icell, jcell, kcell) == -1) THEN
367 imap = imap + 1
368 cellmap(icell, jcell, kcell) = imap
369 cpassert(imap <= nonbonded%nlists)
370 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(imap)
371
372 neighbor_kind_pair%cell_vector(1) = icell
373 neighbor_kind_pair%cell_vector(2) = jcell
374 neighbor_kind_pair%cell_vector(3) = kcell
375 END IF
376 END DO
377 END DO
378 END DO
379 END DO
380 ! Mapping the spherical interaction between subcells
381 ALLOCATE (sphcub(-isubcell(1):isubcell(1), &
382 -isubcell(2):isubcell(2), &
383 -isubcell(3):isubcell(3)))
384 sphcub = .false.
385 IF (all(isubcell /= 0)) THEN
386 radius = real(isubcell(1), kind=dp)**2 + real(isubcell(2), kind=dp)**2 + &
387 REAL(isubcell(3), kind=dp)**2
388 loop1: DO k = -isubcell(3), isubcell(3)
389 loop2: DO j = -isubcell(2), isubcell(2)
390 loop3: DO i = -isubcell(1), isubcell(1)
391 ic = real((/i, j, k/), kind=dp)
392 ! subcell cube vertex
393 DO kx = -1, 1
394 icx(3) = ic(3) + sign(0.5_dp, real(kx, kind=dp))
395 DO jx = -1, 1
396 icx(2) = ic(2) + sign(0.5_dp, real(jx, kind=dp))
397 DO ix = -1, 1
398 icx(1) = ic(1) + sign(0.5_dp, real(ix, kind=dp))
399 vv = icx(1)*icx(1) + icx(2)*icx(2) + icx(3)*icx(3)
400 vv = vv/radius
401 IF (vv <= 1.0_dp) THEN
402 sphcub(i, j, k) = .true.
403 cycle loop3
404 END IF
405 END DO
406 END DO
407 END DO
408 END DO loop3
409 END DO loop2
410 END DO loop1
411 END IF
412 ! Mapping locally all atoms in the zeroth cell
413 ALLOCATE (coord(3, SIZE(particle_set)))
414 DO atom_a = 1, SIZE(particle_set)
415 coord(:, atom_a) = pbc(particle_set(atom_a)%r, cell)
416 END DO
417 ! Associate particles to subcells (local particles)
418 DO ikind = 1, nkind
419 IF (.NOT. ASSOCIATED(atom(ikind)%list_local_a_index)) cycle
420 natom_local_a = SIZE(atom(ikind)%list_local_a_index)
421 DO iatom_local = 1, natom_local_a
422 atom_a = atom(ikind)%list_local_a_index(iatom_local)
423 CALL give_ijk_subcell(coord(:, atom_a), i, j, k, cell, nsubcell)
424 subcell_a(i, j, k)%natom = subcell_a(i, j, k)%natom + 1
425 END DO
426 END DO
427 DO k = 1, nsubcell(3)
428 DO j = 1, nsubcell(2)
429 DO i = 1, nsubcell(1)
430 ALLOCATE (subcell_a(i, j, k)%atom_list(subcell_a(i, j, k)%natom))
431 subcell_a(i, j, k)%natom = 0
432 END DO
433 END DO
434 END DO
435 DO ikind = 1, nkind
436 IF (.NOT. ASSOCIATED(atom(ikind)%list_local_a_index)) cycle
437 natom_local_a = SIZE(atom(ikind)%list_local_a_index)
438 DO iatom_local = 1, natom_local_a
439 atom_a = atom(ikind)%list_local_a_index(iatom_local)
440 CALL give_ijk_subcell(coord(:, atom_a), i, j, k, cell, nsubcell)
441 subcell_a(i, j, k)%natom = subcell_a(i, j, k)%natom + 1
442 subcell_a(i, j, k)%atom_list(subcell_a(i, j, k)%natom) = atom_a
443 END DO
444 END DO
445 ! Associate particles to subcells (distributed particles)
446 DO atom_b = 1, SIZE(particle_set)
447 CALL give_ijk_subcell(coord(:, atom_b), i, j, k, cell, nsubcell)
448 subcell_b(i, j, k)%natom = subcell_b(i, j, k)%natom + 1
449 END DO
450 DO k = 1, nsubcell(3)
451 DO j = 1, nsubcell(2)
452 DO i = 1, nsubcell(1)
453 ALLOCATE (subcell_b(i, j, k)%atom_list(subcell_b(i, j, k)%natom))
454 subcell_b(i, j, k)%natom = 0
455 END DO
456 END DO
457 END DO
458 DO atom_b = 1, SIZE(particle_set)
459 CALL give_ijk_subcell(coord(:, atom_b), i, j, k, cell, nsubcell)
460 subcell_b(i, j, k)%natom = subcell_b(i, j, k)%natom + 1
461 subcell_b(i, j, k)%atom_list(subcell_b(i, j, k)%natom) = atom_b
462 END DO
463 ! Reorder atoms associated to subcells
464 tmpdim = maxval(subcell_a(:, :, :)%natom)
465 tmpdim = max(tmpdim, maxval(subcell_b(:, :, :)%natom))
466 ALLOCATE (work(3*tmpdim))
467 ALLOCATE (kind_of(SIZE(particle_set)))
468 DO i = 1, SIZE(particle_set)
469 kind_of(i) = particle_set(i)%atomic_kind%kind_number
470 END DO
471 DO k = 1, nsubcell(3)
472 DO j = 1, nsubcell(2)
473 DO i = 1, nsubcell(1)
474 CALL reorder_atoms_subcell(subcell_a(i, j, k)%atom_list, kind_of, work)
475 CALL reorder_atoms_subcell(subcell_b(i, j, k)%atom_list, kind_of, work)
476 END DO
477 END DO
478 END DO
479 DEALLOCATE (work, kind_of)
480 zdim = nsubcell(3)
481 ydim = nsubcell(2)
482 xdim = nsubcell(1)
483 is_full = .false.
484 ! We can skip until ik>=0.. this prescreens the order of the subcells
485 ik_start = -isubcell(3)
486 IF (.NOT. any_full) ik_start = 0
487 ! Loop over first subcell
488 loop_a_k: DO a_k = 1, nsubcell(3)
489 loop_a_j: DO a_j = 1, nsubcell(2)
490 loop_a_i: DO a_i = 1, nsubcell(1)
491 IF (subcell_a(a_i, a_j, a_k)%natom == 0) cycle
492 ! Loop over second subcell
493 loop_b_k: DO ik = ik_start, isubcell(3)
494 bg_k = a_k + ik
495 b_k = mod(bg_k, zdim)
496 b_pk = bg_k/zdim
497 IF (b_k <= 0) THEN
498 b_k = zdim + b_k
499 b_pk = b_pk - 1
500 END IF
501 IF ((periodic(3) == 0) .AND. (abs(b_pk) > 0)) cycle
502 ! Setup the starting point.. this prescreens the order of the subcells
503 ij_start = -isubcell(2)
504 IF ((ik == 0) .AND. (ik_start == 0)) ij_start = 0
505 loop_b_j: DO ij = ij_start, isubcell(2)
506 bg_j = a_j + ij
507 b_j = mod(bg_j, ydim)
508 b_pj = bg_j/ydim
509 IF (b_j <= 0) THEN
510 b_j = ydim + b_j
511 b_pj = b_pj - 1
512 END IF
513 IF ((periodic(2) == 0) .AND. (abs(b_pj) > 0)) cycle
514 ! Setup the starting point.. this prescreens the order of the subcells
515 ii_start = -isubcell(1)
516 IF ((ij == 0) .AND. (ij_start == 0)) ii_start = 0
517 loop_b_i: DO ii = ii_start, isubcell(1)
518 ! Ellipsoidal screening of subcells
519 IF (.NOT. sphcub(ii, ij, ik)) cycle
520 bg_i = a_i + ii
521 b_i = mod(bg_i, xdim)
522 b_pi = bg_i/xdim
523 IF (b_i <= 0) THEN
524 b_i = xdim + b_i
525 b_pi = b_pi - 1
526 END IF
527 IF ((periodic(1) == 0) .AND. (abs(b_pi) > 0)) cycle
528 IF (subcell_b(b_i, b_j, b_k)%natom == 0) cycle
529 ! Find the proper neighbor kind pair
530 icellmap = cellmap(b_pi, b_pj, b_pk)
531 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(icellmap)
532 ! Find the replica vector
533 cell_v = 0.0_dp
534 IF ((b_pi /= 0) .OR. (b_pj /= 0) .OR. (b_pk /= 0)) THEN
535 cv_b(1) = b_pi; cv_b(2) = b_pj; cv_b(3) = b_pk
536 CALL scaled_to_real(cell_v, cv_b, cell)
537 END IF
538 subcell000 = (a_k == bg_k) .AND. (a_j == bg_j) .AND. (a_i == bg_i)
539 ! Loop over particles inside subcell_a and subcell_b
540 DO jatom_local = 1, subcell_b(b_i, b_j, b_k)%natom
541 atom_b = subcell_b(b_i, b_j, b_k)%atom_list(jatom_local)
542 jkind = particle_set(atom_b)%atomic_kind%kind_number
543 rb(1) = coord(1, atom_b) + cell_v(1)
544 rb(2) = coord(2, atom_b) + cell_v(2)
545 rb(3) = coord(3, atom_b) + cell_v(3)
546 DO iatom_local = 1, subcell_a(a_i, a_j, a_k)%natom
547 atom_a = subcell_a(a_i, a_j, a_k)%atom_list(iatom_local)
548 ikind = particle_set(atom_a)%atomic_kind%kind_number
549 ! Screen interaction to avoid double counting
550 atom_order = (atom_a <= atom_b)
551 ! Special case for kind combination requiring the full NL
552 IF (any_full) THEN
553 is_full = full_nl(ikind, jkind)
554 IF (is_full) THEN
555 atom_order = (atom_a == atom_b)
556 ELSE
557 IF (ik < 0) cycle
558 IF (ik == 0 .AND. ij < 0) cycle
559 IF (ij == 0 .AND. ii < 0) cycle
560 END IF
561 END IF
562 IF (subcell000 .AND. atom_order) cycle
563 rab(1) = rb(1) - coord(1, atom_a)
564 rab(2) = rb(2) - coord(2, atom_a)
565 rab(3) = rb(3) - coord(3, atom_a)
566 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
567 rab_max = r_max(ikind, jkind)
568 rab2_max = rab_max*rab_max
569 IF (rab2 < rab2_max) THEN
570 ! Diagonal storage
571 j1 = min(ikind, jkind)
572 i1 = max(ikind, jkind) - j1 + 1
573 j1 = nkind - j1 + 1
574 id_kind = nkind00 - (j1*(j1 + 1)/2) + i1
575 ! Store the pair
576 CALL fist_neighbor_add(neighbor_kind_pair, atom_a, atom_b, &
577 rab=rab, &
578 check_spline=check_spline, id_kind=id_kind, &
579 skip=(skip_kind(ikind) .AND. skip_kind(jkind)), &
580 cell=cell, ei_scale14=ei_scale14, &
581 vdw_scale14=vdw_scale14, exclusions=exclusions)
582 ! This is to handle properly when interaction radius is larger than cell size
583 IF ((atom_a == atom_b) .AND. (ik_start == 0)) THEN
584 invcellmap = cellmap(-b_pi, -b_pj, -b_pk)
585 inv_neighbor_kind_pair => nonbonded%neighbor_kind_pairs(invcellmap)
586 rab = rab - 2.0_dp*cell_v
587 CALL fist_neighbor_add(inv_neighbor_kind_pair, atom_a, atom_b, &
588 rab=rab, &
589 check_spline=check_spline, id_kind=id_kind, &
590 skip=(skip_kind(ikind) .AND. skip_kind(jkind)), &
591 cell=cell, ei_scale14=ei_scale14, &
592 vdw_scale14=vdw_scale14, exclusions=exclusions)
593 END IF
594 ! Check for too close hits
595 IF (check_spline) THEN
596 rab2_min = r_minsq(ikind, jkind)
597 IF (rab2 < rab2_min) THEN
599 WRITE (iw, '(T2,A,2I7,2(A,F15.8),A)') "WARNING| Particles: ", &
600 atom_a, atom_b, &
601 " at distance [au]:", sqrt(rab2), " less than: ", &
602 sqrt(rab2_min), &
603 "; increase EMAX_SPLINE."
604 IF (rab2 < rab2_min/(1.06_dp)**2) THEN
605 IF (geo_check) THEN
606 cpabort("GEOMETRY wrong or EMAX_SPLINE too small!")
607 END IF
608 END IF
609 END IF
610 END IF
611 END IF
612 END DO
613 END DO
614 END DO loop_b_i
615 END DO loop_b_j
616 END DO loop_b_k
617 END DO loop_a_i
618 END DO loop_a_j
619 END DO loop_a_k
620 DEALLOCATE (coord)
621 DEALLOCATE (cellmap)
622 DEALLOCATE (sphcub)
623 CALL deallocate_subcell(subcell_a)
624 CALL deallocate_subcell(subcell_b)
625
626 CALL timestop(handle)
627 END SUBROUTINE build_neighbor_lists
628
629! **************************************************************************************************
630!> \brief Write a set of neighbor lists to the output unit.
631!> \param nonbonded ...
632!> \param particle_set ...
633!> \param cell ...
634!> \param para_env ...
635!> \param output_unit ...
636!> \param name ...
637!> \param unit_str ...
638!> \par History
639!> 08.2006 created [tlaino]
640!> \author Teodoro Laino
641! **************************************************************************************************
642 SUBROUTINE write_neighbor_lists(nonbonded, particle_set, cell, para_env, output_unit, &
643 name, unit_str)
644
645 TYPE(fist_neighbor_type), POINTER :: nonbonded
646 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
647 TYPE(cell_type), POINTER :: cell
648 TYPE(mp_para_env_type), POINTER :: para_env
649 INTEGER, INTENT(IN) :: output_unit
650 CHARACTER(LEN=*), INTENT(IN) :: name, unit_str
651
652 CHARACTER(LEN=default_string_length) :: string
653 INTEGER :: atom_a, atom_b, iab, ilist, nneighbor
654 LOGICAL :: print_headline
655 REAL(dp) :: conv, dab
656 REAL(dp), DIMENSION(3) :: cell_v, ra, rab, rb
657 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
658
659 ! Print headline
660 string = ""
661 WRITE (unit=string, fmt="(A,I5,A)") &
662 trim(name)//" IN "//trim(unit_str)//" (PROCESS", para_env%mepos, ")"
663 CALL compress(string)
664 IF (output_unit > 0) WRITE (unit=output_unit, fmt="(/,/,T2,A)") trim(string)
665
666 print_headline = .true.
667 nneighbor = 0
668 conv = cp_unit_from_cp2k(1.0_dp, trim(unit_str))
669 DO iab = 1, SIZE(nonbonded%neighbor_kind_pairs)
670 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(iab)
671 cell_v = matmul(cell%hmat, real(neighbor_kind_pair%cell_vector, kind=dp))
672 DO ilist = 1, neighbor_kind_pair%npairs
673 nneighbor = nneighbor + 1
674 IF (output_unit > 0) THEN
675 ! Print second part of the headline
676 atom_a = neighbor_kind_pair%list(1, ilist)
677 atom_b = neighbor_kind_pair%list(2, ilist)
678 IF (print_headline) THEN
679 WRITE (unit=output_unit, fmt="(T3,2(A6,3(5X,A,5X)),1X,A11,10X,A8,A5,A10,A9)") &
680 "Atom-A", "X", "Y", "Z", "Atom-B", "X", "Y", "Z", "Cell(i,j,k)", &
681 "Distance", "ONFO", "VDW-scale", "EI-scale"
682 print_headline = .false.
683 END IF
684
685 ra(:) = pbc(particle_set(atom_a)%r, cell)
686 rb(:) = pbc(particle_set(atom_b)%r, cell)
687 rab = rb(:) - ra(:) + cell_v
688 dab = sqrt(dot_product(rab, rab))
689 IF (ilist <= neighbor_kind_pair%nscale) THEN
690 WRITE (unit=output_unit, fmt="(T3,2(I6,3(1X,F10.6)),3(1X,I3),10X,F8.4,L4,F11.5,F9.5)") &
691 atom_a, ra(1:3)*conv, &
692 atom_b, rb(1:3)*conv, &
693 neighbor_kind_pair%cell_vector, &
694 dab*conv, &
695 neighbor_kind_pair%is_onfo(ilist), &
696 neighbor_kind_pair%vdw_scale(ilist), &
697 neighbor_kind_pair%ei_scale(ilist)
698 ELSE
699 WRITE (unit=output_unit, fmt="(T3,2(I6,3(1X,F10.6)),3(1X,I3),10X,F8.4)") &
700 atom_a, ra(1:3)*conv, &
701 atom_b, rb(1:3)*conv, &
702 neighbor_kind_pair%cell_vector, &
703 dab*conv
704 END IF
705 END IF
706 END DO ! ilist
707 END DO ! iab
708
709 string = ""
710 WRITE (unit=string, fmt="(A,I12,A,I12)") &
711 "Total number of neighbor interactions for process", para_env%mepos, ":", &
712 nneighbor
713 CALL compress(string)
714 IF (output_unit > 0) WRITE (unit=output_unit, fmt="(/,T2,A)") trim(string)
715
716 END SUBROUTINE write_neighbor_lists
717
718! **************************************************************************************************
719!> \brief Sort the generated neighbor list according the kind
720!> \param nonbonded ...
721!> \param nkinds ...
722!> \par History
723!> 09.2007 created [tlaino] University of Zurich - Reducing memory usage
724!> for the FIST neighbor lists
725!> \author Teodoro Laino - University of Zurich
726! **************************************************************************************************
727 SUBROUTINE sort_neighbor_lists(nonbonded, nkinds)
728
729 TYPE(fist_neighbor_type), POINTER :: nonbonded
730 INTEGER, INTENT(IN) :: nkinds
731
732 CHARACTER(LEN=*), PARAMETER :: routinen = 'sort_neighbor_lists'
733
734 INTEGER :: handle, iab, id_kind, ikind, ipair, &
735 jkind, max_alloc_size, npairs, nscale, &
736 tmp
737 INTEGER, ALLOCATABLE, DIMENSION(:) :: indj
738 INTEGER, DIMENSION(:), POINTER :: work
739 INTEGER, DIMENSION(:, :), POINTER :: list_copy
740 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
741
742 NULLIFY (neighbor_kind_pair)
743 CALL timeset(routinen, handle)
744 ! define a lookup table to get jkind for a given id_kind
745 ALLOCATE (indj(nkinds*(nkinds + 1)/2))
746 id_kind = 0
747 DO jkind = 1, nkinds
748 DO ikind = jkind, nkinds
749 id_kind = id_kind + 1
750 indj(id_kind) = jkind
751 END DO
752 END DO
753 ! loop over all nlists and sort the pairs within each list.
754 DO iab = 1, nonbonded%nlists
755 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(iab)
756 npairs = neighbor_kind_pair%npairs
757 nscale = neighbor_kind_pair%nscale
758 IF (npairs /= 0) THEN
759 IF (npairs > nscale) THEN
760 ! 1) Sort the atom pairs by id_kind. Pairs whose interactions are
761 ! scaled (possibly to zero for exclusion) are not touched. They
762 ! stay packed in the beginning. Sorting is skipped altogether when
763 ! all pairs have scaled interactions.
764 ALLOCATE (work(1:npairs - nscale))
765 ALLOCATE (list_copy(2, 1:npairs - nscale))
766 ! Copy of the pair list is required to perform the permutation below
767 ! correctly.
768 list_copy = neighbor_kind_pair%list(:, nscale + 1:npairs)
769 CALL sort(neighbor_kind_pair%id_kind(nscale + 1:npairs), npairs - nscale, work)
770 ! Reorder atoms using the same permutation that was used to sort
771 ! the array id_kind.
772 DO ipair = nscale + 1, npairs
773 tmp = work(ipair - nscale)
774 neighbor_kind_pair%list(1, ipair) = list_copy(1, tmp)
775 neighbor_kind_pair%list(2, ipair) = list_copy(2, tmp)
776 END DO
777 DEALLOCATE (work)
778 DEALLOCATE (list_copy)
779 END IF
780 ! 2) determine the intervals (groups) in the pair list that correspond
781 ! to a certain id_kind. also store the corresponding ikind and
782 ! jkind. Note that this part does not assume ikind to be sorted,
783 ! but it only makes sense when contiguous blobs of the same ikind
784 ! are present.
785 ! Allocate sufficient memory in case all pairs of atom kinds are
786 ! present, and also provide storage for the pairs with exclusion
787 ! flags, which are unsorted.
788 max_alloc_size = nkinds*(nkinds + 1)/2 + nscale
789 IF (ASSOCIATED(neighbor_kind_pair%grp_kind_start)) THEN
790 DEALLOCATE (neighbor_kind_pair%grp_kind_start)
791 END IF
792 ALLOCATE (neighbor_kind_pair%grp_kind_start(max_alloc_size))
793 IF (ASSOCIATED(neighbor_kind_pair%grp_kind_end)) THEN
794 DEALLOCATE (neighbor_kind_pair%grp_kind_end)
795 END IF
796 ALLOCATE (neighbor_kind_pair%grp_kind_end(max_alloc_size))
797 IF (ASSOCIATED(neighbor_kind_pair%ij_kind)) THEN
798 DEALLOCATE (neighbor_kind_pair%ij_kind)
799 END IF
800 ALLOCATE (neighbor_kind_pair%ij_kind(2, max_alloc_size))
801 ! Start the first interval.
802 ipair = 1
803 neighbor_kind_pair%ngrp_kind = 1
804 neighbor_kind_pair%grp_kind_start(neighbor_kind_pair%ngrp_kind) = ipair
805 ! Get ikind and jkind corresponding to id_kind.
806 id_kind = neighbor_kind_pair%id_kind(ipair)
807 jkind = indj(id_kind)
808 tmp = nkinds - jkind
809 ikind = nkinds + id_kind - nkinds*(nkinds + 1)/2 + (tmp*(tmp + 1)/2)
810 neighbor_kind_pair%ij_kind(1, neighbor_kind_pair%ngrp_kind) = ikind
811 neighbor_kind_pair%ij_kind(2, neighbor_kind_pair%ngrp_kind) = jkind
812 ! Define the remaining intervals.
813 DO ipair = 2, npairs
814 IF (neighbor_kind_pair%id_kind(ipair) /= neighbor_kind_pair%id_kind(ipair - 1)) THEN
815 neighbor_kind_pair%grp_kind_end(neighbor_kind_pair%ngrp_kind) = ipair - 1
816 neighbor_kind_pair%ngrp_kind = neighbor_kind_pair%ngrp_kind + 1
817 neighbor_kind_pair%grp_kind_start(neighbor_kind_pair%ngrp_kind) = ipair
818 ! Get ikind and jkind corresponding to id_kind.
819 id_kind = neighbor_kind_pair%id_kind(ipair)
820 jkind = indj(id_kind)
821 tmp = nkinds - jkind
822 ikind = nkinds + id_kind - nkinds*(nkinds + 1)/2 + (tmp*(tmp + 1)/2)
823 neighbor_kind_pair%ij_kind(1, neighbor_kind_pair%ngrp_kind) = ikind
824 neighbor_kind_pair%ij_kind(2, neighbor_kind_pair%ngrp_kind) = jkind
825 END IF
826 END DO
827 ! Finish the last interval.
828 neighbor_kind_pair%grp_kind_end(neighbor_kind_pair%ngrp_kind) = npairs
829 ! Reduce the grp arrays to the actual size because not all pairs of
830 ! atom types have to be present in this pair list.
831 CALL reallocate(neighbor_kind_pair%grp_kind_start, 1, neighbor_kind_pair%ngrp_kind)
832 CALL reallocate(neighbor_kind_pair%grp_kind_end, 1, neighbor_kind_pair%ngrp_kind)
833 CALL reallocate(neighbor_kind_pair%ij_kind, 1, 2, 1, neighbor_kind_pair%ngrp_kind)
834 END IF
835 ! Clean the memory..
836 DEALLOCATE (neighbor_kind_pair%id_kind)
837 END DO
838 DEALLOCATE (indj)
839 CALL timestop(handle)
840 END SUBROUTINE sort_neighbor_lists
841
842END MODULE fist_neighbor_lists
Definition atom.F:9
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Handles all functions related to the CELL.
Definition cell_types.F:15
subroutine, public scaled_to_real(r, s, cell)
Transform scaled cell coordinates real coordinates. r=h*s.
Definition cell_types.F:516
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:195
real(kind=dp) function, public plane_distance(h, k, l, cell)
Calculate the distance between two lattice planes as defined by a triple of Miller indices (hkl).
Definition cell_types.F:252
various routines to log and control the output. The idea is that decisions about where to log should ...
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
unit conversion facility
Definition cp_units.F:30
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1179
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
an exclusion type
Define the neighbor list data types and the corresponding functionality.
subroutine, public fist_neighbor_add(neighbor_kind_pair, atom_a, atom_b, rab, check_spline, id_kind, skip, cell, ei_scale14, vdw_scale14, exclusions)
...
subroutine, public fist_neighbor_init(fist_neighbor, ncell)
...
subroutine, public fist_neighbor_deallocate(fist_neighbor)
...
Generate the atomic neighbor lists for FIST.
subroutine, public build_fist_neighbor_lists(atomic_kind_set, particle_set, local_particles, cell, r_max, r_minsq, ei_scale14, vdw_scale14, nonbonded, para_env, build_from_scratch, geo_check, mm_section, full_nl, exclusions)
...
objects that represent the structure of input sections and the data contained in an input section
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
Utility routines for the memory handling.
Interface to the message passing library MPI.
Define the data structure for the particle information.
logical function, public qmmm_ff_precond_only_qm(id1, id2, id3, id4, is_link)
This function handles the atom names and modifies the "_QM_" prefix, in order to find the parameters ...
Utilities for string manipulations.
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.
subcell types and allocation routines
subroutine, public deallocate_subcell(subcell)
Deallocate a subcell grid structure.
subroutine, public give_ijk_subcell(r, i, j, k, cell, nsubcell)
...
subroutine, public reorder_atoms_subcell(atom_list, kind_of, work)
...
subroutine, public allocate_subcell(subcell, nsubcell, maxatom, cell)
Allocate and initialize a subcell grid structure for the atomic neighbor search.
All kind of helpful little routines.
Definition util.F:14
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
type of a logger, at the moment it contains just a print level starting at which level it should be l...
structure to store local (to a processor) ordered lists of integers.
A type used to store lists of exclusions and onfos.
stores all the informations relevant to an mpi environment