(git:1f285aa)
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 
21  USE atomic_kind_types, ONLY: atomic_kind_type,&
24  USE cell_types, ONLY: cell_type,&
25  get_cell,&
26  pbc,&
31  cp_logger_type
32  USE cp_output_handling, ONLY: cp_p_file,&
36  USE cp_units, ONLY: cp_unit_from_cp2k
37  USE distribution_1d_types, ONLY: distribution_1d_type
38  USE exclusion_types, ONLY: exclusion_type
42  fist_neighbor_type,&
43  neighbor_kind_pairs_type
44  USE input_section_types, ONLY: section_vals_type,&
46  USE kinds, ONLY: default_string_length,&
47  dp
48  USE memory_utilities, ONLY: reallocate
49  USE message_passing, ONLY: mp_para_env_type
50  USE particle_types, ONLY: particle_type
52  USE string_utilities, ONLY: compress
53  USE subcell_types, ONLY: allocate_subcell,&
57  subcell_type
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 
76 CONTAINS
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 
842 END MODULE fist_neighbor_lists
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
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 ...
Definition: qmmm_ff_fist.F:39
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
Definition: subcell_types.F:14
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.
Definition: subcell_types.F:55
All kind of helpful little routines.
Definition: util.F:14