Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 : ! **************************************************************************************************
19 : MODULE fist_neighbor_lists
20 :
21 : USE atomic_kind_types, ONLY: atomic_kind_type,&
22 : get_atomic_kind,&
23 : get_atomic_kind_set
24 : USE cell_types, ONLY: cell_type,&
25 : get_cell,&
26 : pbc,&
27 : plane_distance,&
28 : scaled_to_real
29 : USE cp_log_handling, ONLY: cp_get_default_logger,&
30 : cp_logger_get_default_unit_nr,&
31 : cp_logger_type
32 : USE cp_output_handling, ONLY: cp_p_file,&
33 : cp_print_key_finished_output,&
34 : cp_print_key_should_output,&
35 : cp_print_key_unit_nr
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
39 : USE fist_neighbor_list_types, ONLY: fist_neighbor_add,&
40 : fist_neighbor_deallocate,&
41 : fist_neighbor_init,&
42 : fist_neighbor_type,&
43 : neighbor_kind_pairs_type
44 : USE input_section_types, ONLY: section_vals_type,&
45 : section_vals_val_get
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
51 : USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm
52 : USE string_utilities, ONLY: compress
53 : USE subcell_types, ONLY: allocate_subcell,&
54 : deallocate_subcell,&
55 : give_ijk_subcell,&
56 : reorder_atoms_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 => NULL(), &
70 : list_local_a_index => NULL()
71 : END TYPE local_atoms_type
72 :
73 : ! Public subroutines
74 : PUBLIC :: build_fist_neighbor_lists
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 19326 : SUBROUTINE build_fist_neighbor_lists(atomic_kind_set, particle_set, &
100 19326 : local_particles, cell, r_max, r_minsq, ei_scale14, vdw_scale14, &
101 : nonbonded, para_env, build_from_scratch, geo_check, mm_section, &
102 19326 : 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 19326 : LOGICAL, DIMENSION(:), POINTER :: skip_kind
126 19326 : 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 19326 : CALL timeset(routineN, handle)
132 19326 : NULLIFY (logger)
133 19326 : logger => cp_get_default_logger()
134 :
135 19326 : print_subcell_grid = .FALSE.
136 : output_unit = cp_print_key_unit_nr(logger, mm_section, "PRINT%SUBCELL", &
137 19326 : extension=".Log")
138 19326 : IF (output_unit > 0) print_subcell_grid = .TRUE.
139 :
140 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
141 19326 : maxatom=maxatom)
142 :
143 19326 : 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 19326 : IF (PRESENT(exclusions)) THEN
148 10839 : CPASSERT(present_local_particles)
149 : END IF
150 :
151 : ! Allocate work storage
152 19326 : nkind = SIZE(atomic_kind_set)
153 118257 : ALLOCATE (atom(nkind))
154 57978 : ALLOCATE (skip_kind(nkind))
155 : ! full_nl
156 19326 : IF (PRESENT(full_nl)) THEN
157 11171 : my_full_nl => full_nl
158 : ELSE
159 32620 : ALLOCATE (my_full_nl(nkind, nkind))
160 496909 : my_full_nl = .FALSE.
161 : END IF
162 : ! Initialize the local data structures
163 79605 : DO ikind = 1, nkind
164 60279 : atomic_kind => atomic_kind_set(ikind)
165 60279 : NULLIFY (atom(ikind)%list)
166 60279 : NULLIFY (atom(ikind)%list_local_a_index)
167 :
168 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
169 60279 : atom_list=atom(ikind)%list, name=kind_name)
170 60279 : skip_kind(ikind) = qmmm_ff_precond_only_qm(kind_name)
171 60279 : IF (present_local_particles) THEN
172 50224 : natom_local_a = local_particles%n_el(ikind)
173 : ELSE
174 10055 : natom_local_a = SIZE(atom(ikind)%list)
175 : END IF
176 79605 : IF (natom_local_a > 0) THEN
177 167247 : ALLOCATE (atom(ikind)%list_local_a_index(natom_local_a))
178 : ! Build index vector for mapping
179 1075400 : DO iatom_local = 1, natom_local_a
180 1019651 : IF (present_local_particles) THEN
181 703427 : atom_a = local_particles%list(ikind)%array(iatom_local)
182 : ELSE
183 316224 : atom_a = atom(ikind)%list(iatom_local)
184 : END IF
185 1075400 : atom(ikind)%list_local_a_index(iatom_local) = atom_a
186 : END DO
187 :
188 : END IF
189 : END DO
190 :
191 19326 : IF (build_from_scratch) THEN
192 10660 : IF (ASSOCIATED(nonbonded)) THEN
193 0 : 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 27813 : my_full_nl, exclusions)
202 :
203 : ! Sort the list according kinds for each cell
204 19326 : CALL sort_neighbor_lists(nonbonded, nkind)
205 :
206 19326 : print_key_path = "PRINT%NEIGHBOR_LISTS"
207 :
208 19326 : 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 460 : file_position="REWIND")
218 460 : 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 460 : "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 460 : local=.TRUE.)
226 : END IF
227 :
228 : ! Release work storage
229 79605 : DO ikind = 1, nkind
230 60279 : NULLIFY (atom(ikind)%list)
231 79605 : IF (ASSOCIATED(atom(ikind)%list_local_a_index)) THEN
232 55749 : DEALLOCATE (atom(ikind)%list_local_a_index)
233 : END IF
234 : END DO
235 19326 : IF (PRESENT(full_nl)) THEN
236 : NULLIFY (my_full_nl)
237 : ELSE
238 8155 : DEALLOCATE (my_full_nl)
239 : END IF
240 19326 : DEALLOCATE (atom)
241 :
242 19326 : DEALLOCATE (skip_kind)
243 :
244 : CALL cp_print_key_finished_output(unit_nr=output_unit, &
245 : logger=logger, &
246 : basis_section=mm_section, &
247 19326 : print_key_path="PRINT%SUBCELL")
248 19326 : CALL timestop(handle)
249 :
250 38652 : 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 19326 : SUBROUTINE build_neighbor_lists(nonbonded, particle_set, atom, cell, &
274 38652 : print_subcell_grid, output_unit, r_max, r_minsq, &
275 19326 : 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 19326 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, work
298 19326 : 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 19326 : 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 19326 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coord
307 : TYPE(neighbor_kind_pairs_type), POINTER :: inv_neighbor_kind_pair, &
308 : neighbor_kind_pair
309 19326 : TYPE(subcell_type), DIMENSION(:, :, :), POINTER :: subcell_a, subcell_b
310 :
311 19326 : CALL timeset(routineN, handle)
312 19326 : nkind = SIZE(atom)
313 77304 : nsubcell = 1
314 19326 : isubcell = 0
315 19326 : ncell = 0
316 1865161 : any_full = ANY(full_nl)
317 : CALL get_cell(cell=cell, &
318 : abc=abc, &
319 19326 : periodic=periodic)
320 : ! Determines the number of subcells
321 79605 : DO ikind = 1, nkind
322 1004169 : DO jkind = ikind, nkind
323 : ! Calculate the square of the maximum interaction distance
324 924564 : rab_max = r_max(ikind, jkind)
325 924564 : IF (skip_kind(ikind) .AND. skip_kind(jkind)) CYCLE
326 912106 : nsubcell(1) = MAX(nsubcell(1), CEILING(plane_distance(1, 0, 0, cell)/rab_max))
327 912106 : nsubcell(2) = MAX(nsubcell(2), CEILING(plane_distance(0, 1, 0, cell)/rab_max))
328 984843 : 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 79605 : DO ikind = 1, nkind
333 1004169 : DO jkind = ikind, nkind
334 924564 : IF (skip_kind(ikind) .AND. skip_kind(jkind)) CYCLE
335 : ! Calculate the square of the maximum interaction distance
336 912106 : rab_max = r_max(ikind, jkind)
337 912106 : sab_max(1) = rab_max/plane_distance(1, 0, 0, cell)
338 912106 : sab_max(2) = rab_max/plane_distance(0, 1, 0, cell)
339 912106 : sab_max(3) = rab_max/plane_distance(0, 0, 1, cell)
340 3648424 : ncell = MAX(ncell(:), CEILING(sab_max(:)*periodic(:)))
341 3721161 : isubcell = MAX(isubcell(:), CEILING(sab_max(:)*REAL(nsubcell(:), KIND=dp)))
342 : END DO
343 : END DO
344 19326 : CALL fist_neighbor_init(nonbonded, ncell)
345 : ! Print headline
346 19326 : IF (print_subcell_grid) THEN
347 : WRITE (UNIT=output_unit, FMT="(/,/,T2,A,/)") &
348 86 : "SUBCELL GRID INFO FOR THE "//TRIM(name)//" NEIGHBOR LISTS"
349 86 : WRITE (UNIT=output_unit, FMT="(T4,A,10X,3I10)") " NUMBER OF SUBCELLS ::", nsubcell
350 86 : WRITE (UNIT=output_unit, FMT="(T4,A,10X,3I10)") " NUMBER OF PERIODIC IMAGES ::", ncell
351 86 : WRITE (UNIT=output_unit, FMT="(T4,A,10X,3I10)") " NUMBER OF INTERACTING SUBCELLS ::", isubcell
352 : END IF
353 : ! Allocate subcells
354 19326 : CALL allocate_subcell(subcell_a, nsubcell, cell=cell)
355 19326 : CALL allocate_subcell(subcell_b, nsubcell, cell=cell)
356 : ! Let's map the sequence of the periodic images
357 77304 : ncellmax = MAXVAL(ncell)
358 96630 : ALLOCATE (cellmap(-ncellmax:ncellmax, -ncellmax:ncellmax, -ncellmax:ncellmax))
359 1115660 : cellmap = -1
360 19326 : imap = 0
361 19326 : nkind00 = nkind*(nkind + 1)/2
362 56197 : DO imax_cell = 0, ncellmax
363 134076 : DO kcell = -imax_cell, imax_cell
364 367661 : DO jcell = -imax_cell, imax_cell
365 1867845 : DO icell = -imax_cell, imax_cell
366 1789966 : IF (cellmap(icell, jcell, kcell) == -1) THEN
367 858560 : imap = imap + 1
368 858560 : cellmap(icell, jcell, kcell) = imap
369 858560 : CPASSERT(imap <= nonbonded%nlists)
370 858560 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(imap)
371 :
372 858560 : neighbor_kind_pair%cell_vector(1) = icell
373 858560 : neighbor_kind_pair%cell_vector(2) = jcell
374 858560 : 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 96630 : -isubcell(3):isubcell(3)))
384 3113004 : sphcub = .FALSE.
385 77298 : IF (ALL(isubcell /= 0)) THEN
386 : radius = REAL(isubcell(1), KIND=dp)**2 + REAL(isubcell(2), KIND=dp)**2 + &
387 19324 : REAL(isubcell(3), KIND=dp)**2
388 113340 : loop1: DO k = -isubcell(3), isubcell(3)
389 584628 : loop2: DO j = -isubcell(2), isubcell(2)
390 3093672 : loop3: DO i = -isubcell(1), isubcell(1)
391 10113472 : ic = REAL((/i, j, k/), KIND=dp)
392 : ! subcell cube vertex
393 3000034 : DO kx = -1, 1
394 2528746 : icx(3) = ic(3) + SIGN(0.5_dp, REAL(kx, KIND=dp))
395 2622052 : DO jx = -1, 1
396 2621674 : icx(2) = ic(2) + SIGN(0.5_dp, REAL(jx, KIND=dp))
397 3068784 : DO ix = -1, 1
398 2975100 : icx(1) = ic(1) + SIGN(0.5_dp, REAL(ix, KIND=dp))
399 2975100 : vv = icx(1)*icx(1) + icx(2)*icx(2) + icx(3)*icx(3)
400 2975100 : vv = vv/radius
401 3068406 : IF (vv <= 1.0_dp) THEN
402 2528368 : sphcub(i, j, k) = .TRUE.
403 2528368 : 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 57978 : ALLOCATE (coord(3, SIZE(particle_set)))
414 1501066 : DO atom_a = 1, SIZE(particle_set)
415 1501066 : coord(:, atom_a) = pbc(particle_set(atom_a)%r, cell)
416 : END DO
417 : ! Associate particles to subcells (local particles)
418 79605 : DO ikind = 1, nkind
419 60279 : IF (.NOT. ASSOCIATED(atom(ikind)%list_local_a_index)) CYCLE
420 55749 : natom_local_a = SIZE(atom(ikind)%list_local_a_index)
421 1094726 : DO iatom_local = 1, natom_local_a
422 1019651 : atom_a = atom(ikind)%list_local_a_index(iatom_local)
423 1019651 : CALL give_ijk_subcell(coord(:, atom_a), i, j, k, cell, nsubcell)
424 1079930 : subcell_a(i, j, k)%natom = subcell_a(i, j, k)%natom + 1
425 : END DO
426 : END DO
427 78480 : DO k = 1, nsubcell(3)
428 336705 : DO j = 1, nsubcell(2)
429 2272858 : DO i = 1, nsubcell(1)
430 4099529 : ALLOCATE (subcell_a(i, j, k)%atom_list(subcell_a(i, j, k)%natom))
431 2213704 : subcell_a(i, j, k)%natom = 0
432 : END DO
433 : END DO
434 : END DO
435 79605 : DO ikind = 1, nkind
436 60279 : IF (.NOT. ASSOCIATED(atom(ikind)%list_local_a_index)) CYCLE
437 55749 : natom_local_a = SIZE(atom(ikind)%list_local_a_index)
438 1094726 : DO iatom_local = 1, natom_local_a
439 1019651 : atom_a = atom(ikind)%list_local_a_index(iatom_local)
440 1019651 : CALL give_ijk_subcell(coord(:, atom_a), i, j, k, cell, nsubcell)
441 1019651 : subcell_a(i, j, k)%natom = subcell_a(i, j, k)%natom + 1
442 1079930 : 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 1501066 : DO atom_b = 1, SIZE(particle_set)
447 1481740 : CALL give_ijk_subcell(coord(:, atom_b), i, j, k, cell, nsubcell)
448 1501066 : subcell_b(i, j, k)%natom = subcell_b(i, j, k)%natom + 1
449 : END DO
450 78480 : DO k = 1, nsubcell(3)
451 336705 : DO j = 1, nsubcell(2)
452 2272858 : DO i = 1, nsubcell(1)
453 4105981 : ALLOCATE (subcell_b(i, j, k)%atom_list(subcell_b(i, j, k)%natom))
454 2213704 : subcell_b(i, j, k)%natom = 0
455 : END DO
456 : END DO
457 : END DO
458 1501066 : DO atom_b = 1, SIZE(particle_set)
459 1481740 : CALL give_ijk_subcell(coord(:, atom_b), i, j, k, cell, nsubcell)
460 1481740 : subcell_b(i, j, k)%natom = subcell_b(i, j, k)%natom + 1
461 1501066 : 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 2292184 : tmpdim = MAXVAL(subcell_a(:, :, :)%natom)
465 2292184 : tmpdim = MAX(tmpdim, MAXVAL(subcell_b(:, :, :)%natom))
466 57978 : ALLOCATE (work(3*tmpdim))
467 57978 : ALLOCATE (kind_of(SIZE(particle_set)))
468 1501066 : DO i = 1, SIZE(particle_set)
469 1501066 : kind_of(i) = particle_set(i)%atomic_kind%kind_number
470 : END DO
471 78480 : DO k = 1, nsubcell(3)
472 336705 : DO j = 1, nsubcell(2)
473 2272858 : DO i = 1, nsubcell(1)
474 1955479 : CALL reorder_atoms_subcell(subcell_a(i, j, k)%atom_list, kind_of, work)
475 2213704 : CALL reorder_atoms_subcell(subcell_b(i, j, k)%atom_list, kind_of, work)
476 : END DO
477 : END DO
478 : END DO
479 19326 : DEALLOCATE (work, kind_of)
480 19326 : zdim = nsubcell(3)
481 19326 : ydim = nsubcell(2)
482 19326 : xdim = nsubcell(1)
483 19326 : is_full = .FALSE.
484 : ! We can skip until ik>=0.. this prescreens the order of the subcells
485 19326 : ik_start = -isubcell(3)
486 19326 : IF (.NOT. any_full) ik_start = 0
487 : ! Loop over first subcell
488 78480 : loop_a_k: DO a_k = 1, nsubcell(3)
489 336705 : loop_a_j: DO a_j = 1, nsubcell(2)
490 2272858 : loop_a_i: DO a_i = 1, nsubcell(1)
491 1955479 : IF (subcell_a(a_i, a_j, a_k)%natom == 0) CYCLE
492 : ! Loop over second subcell
493 1008890 : loop_b_k: DO ik = ik_start, isubcell(3)
494 562094 : bg_k = a_k + ik
495 562094 : b_k = MOD(bg_k, zdim)
496 562094 : b_pk = bg_k/zdim
497 562094 : IF (b_k <= 0) THEN
498 123788 : b_k = zdim + b_k
499 123788 : b_pk = b_pk - 1
500 : END IF
501 562094 : IF ((periodic(3) == 0) .AND. (ABS(b_pk) > 0)) CYCLE
502 : ! Setup the starting point.. this prescreens the order of the subcells
503 554031 : ij_start = -isubcell(2)
504 554031 : IF ((ik == 0) .AND. (ik_start == 0)) ij_start = 0
505 4893012 : loop_b_j: DO ij = ij_start, isubcell(2)
506 2383502 : bg_j = a_j + ij
507 2383502 : b_j = MOD(bg_j, ydim)
508 2383502 : b_pj = bg_j/ydim
509 2383502 : IF (b_j <= 0) THEN
510 600156 : b_j = ydim + b_j
511 600156 : b_pj = b_pj - 1
512 : END IF
513 2383502 : IF ((periodic(2) == 0) .AND. (ABS(b_pj) > 0)) CYCLE
514 : ! Setup the starting point.. this prescreens the order of the subcells
515 2367348 : ii_start = -isubcell(1)
516 2367348 : IF ((ij == 0) .AND. (ij_start == 0)) ii_start = 0
517 14438039 : loop_b_i: DO ii = ii_start, isubcell(1)
518 : ! Ellipsoidal screening of subcells
519 11508597 : IF (.NOT. sphcub(ii, ij, ik)) CYCLE
520 11508596 : bg_i = a_i + ii
521 11508596 : b_i = MOD(bg_i, xdim)
522 11508596 : b_pi = bg_i/xdim
523 11508596 : IF (b_i <= 0) THEN
524 3077687 : b_i = xdim + b_i
525 3077687 : b_pi = b_pi - 1
526 : END IF
527 11508596 : IF ((periodic(1) == 0) .AND. (ABS(b_pi) > 0)) CYCLE
528 11470469 : IF (subcell_b(b_i, b_j, b_k)%natom == 0) CYCLE
529 : ! Find the proper neighbor kind pair
530 8106774 : icellmap = cellmap(b_pi, b_pj, b_pk)
531 8106774 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(icellmap)
532 : ! Find the replica vector
533 8106774 : cell_v = 0.0_dp
534 8106774 : IF ((b_pi /= 0) .OR. (b_pj /= 0) .OR. (b_pk /= 0)) THEN
535 3854334 : cv_b(1) = b_pi; cv_b(2) = b_pj; cv_b(3) = b_pk
536 3854334 : CALL scaled_to_real(cell_v, cv_b, cell)
537 : END IF
538 8106774 : 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 93317802 : DO jatom_local = 1, subcell_b(b_i, b_j, b_k)%natom
541 82827526 : atom_b = subcell_b(b_i, b_j, b_k)%atom_list(jatom_local)
542 82827526 : jkind = particle_set(atom_b)%atomic_kind%kind_number
543 82827526 : rb(1) = coord(1, atom_b) + cell_v(1)
544 82827526 : rb(2) = coord(2, atom_b) + cell_v(2)
545 82827526 : rb(3) = coord(3, atom_b) + cell_v(3)
546 2488256709 : DO iatom_local = 1, subcell_a(a_i, a_j, a_k)%natom
547 2393920586 : atom_a = subcell_a(a_i, a_j, a_k)%atom_list(iatom_local)
548 2393920586 : ikind = particle_set(atom_a)%atomic_kind%kind_number
549 : ! Screen interaction to avoid double counting
550 2393920586 : atom_order = (atom_a <= atom_b)
551 : ! Special case for kind combination requiring the full NL
552 2393920586 : IF (any_full) THEN
553 19336474 : is_full = full_nl(ikind, jkind)
554 19336474 : IF (is_full) THEN
555 9066349 : atom_order = (atom_a == atom_b)
556 : ELSE
557 10270125 : IF (ik < 0) CYCLE
558 6900294 : IF (ik == 0 .AND. ij < 0) CYCLE
559 5491910 : IF (ij == 0 .AND. ii < 0) CYCLE
560 : END IF
561 : END IF
562 2388586438 : IF (subcell000 .AND. atom_order) CYCLE
563 2368390332 : rab(1) = rb(1) - coord(1, atom_a)
564 2368390332 : rab(2) = rb(2) - coord(2, atom_a)
565 2368390332 : rab(3) = rb(3) - coord(3, atom_a)
566 2368390332 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
567 2368390332 : rab_max = r_max(ikind, jkind)
568 2368390332 : rab2_max = rab_max*rab_max
569 2451217858 : IF (rab2 < rab2_max) THEN
570 : ! Diagonal storage
571 151495343 : j1 = MIN(ikind, jkind)
572 151495343 : i1 = MAX(ikind, jkind) - j1 + 1
573 151495343 : j1 = nkind - j1 + 1
574 151495343 : 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 310480404 : vdw_scale14=vdw_scale14, exclusions=exclusions)
582 : ! This is to handle properly when interaction radius is larger than cell size
583 151495343 : IF ((atom_a == atom_b) .AND. (ik_start == 0)) THEN
584 352913 : invcellmap = cellmap(-b_pi, -b_pj, -b_pk)
585 352913 : inv_neighbor_kind_pair => nonbonded%neighbor_kind_pairs(invcellmap)
586 1411652 : 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 725693 : vdw_scale14=vdw_scale14, exclusions=exclusions)
593 : END IF
594 : ! Check for too close hits
595 151495343 : IF (check_spline) THEN
596 150834642 : rab2_min = r_minsq(ikind, jkind)
597 150834642 : IF (rab2 < rab2_min) THEN
598 1 : iw = cp_logger_get_default_unit_nr()
599 1 : WRITE (iw, '(T2,A,2I7,2(A,F15.8),A)') "WARNING| Particles: ", &
600 1 : atom_a, atom_b, &
601 1 : " at distance [au]:", SQRT(rab2), " less than: ", &
602 1 : SQRT(rab2_min), &
603 2 : "; increase EMAX_SPLINE."
604 1 : IF (rab2 < rab2_min/(1.06_dp)**2) THEN
605 0 : IF (geo_check) THEN
606 0 : 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 19326 : DEALLOCATE (coord)
621 19326 : DEALLOCATE (cellmap)
622 19326 : DEALLOCATE (sphcub)
623 19326 : CALL deallocate_subcell(subcell_a)
624 19326 : CALL deallocate_subcell(subcell_b)
625 :
626 19326 : CALL timestop(handle)
627 19326 : 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 460 : 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 460 : string = ""
661 : WRITE (UNIT=string, FMT="(A,I5,A)") &
662 460 : TRIM(name)//" IN "//TRIM(unit_str)//" (PROCESS", para_env%mepos, ")"
663 460 : CALL compress(string)
664 460 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") TRIM(string)
665 :
666 460 : print_headline = .TRUE.
667 460 : nneighbor = 0
668 460 : conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
669 22116 : DO iab = 1, SIZE(nonbonded%neighbor_kind_pairs)
670 21656 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(iab)
671 346496 : cell_v = MATMUL(cell%hmat, REAL(neighbor_kind_pair%cell_vector, KIND=dp))
672 274699 : DO ilist = 1, neighbor_kind_pair%npairs
673 252583 : nneighbor = nneighbor + 1
674 274239 : IF (output_unit > 0) THEN
675 : ! Print second part of the headline
676 252583 : atom_a = neighbor_kind_pair%list(1, ilist)
677 252583 : atom_b = neighbor_kind_pair%list(2, ilist)
678 252583 : 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 294 : "Atom-A", "X", "Y", "Z", "Atom-B", "X", "Y", "Z", "Cell(i,j,k)", &
681 588 : "Distance", "ONFO", "VDW-scale", "EI-scale"
682 294 : print_headline = .FALSE.
683 : END IF
684 :
685 252583 : ra(:) = pbc(particle_set(atom_a)%r, cell)
686 252583 : rb(:) = pbc(particle_set(atom_b)%r, cell)
687 1010332 : rab = rb(:) - ra(:) + cell_v
688 1010332 : dab = SQRT(DOT_PRODUCT(rab, rab))
689 252583 : 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 27248 : atom_a, ra(1:3)*conv, &
692 27248 : atom_b, rb(1:3)*conv, &
693 6812 : neighbor_kind_pair%cell_vector, &
694 6812 : dab*conv, &
695 6812 : neighbor_kind_pair%is_onfo(ilist), &
696 6812 : neighbor_kind_pair%vdw_scale(ilist), &
697 13624 : 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 983084 : atom_a, ra(1:3)*conv, &
701 983084 : atom_b, rb(1:3)*conv, &
702 245771 : neighbor_kind_pair%cell_vector, &
703 491542 : dab*conv
704 : END IF
705 : END IF
706 : END DO ! ilist
707 : END DO ! iab
708 :
709 460 : string = ""
710 : WRITE (UNIT=string, FMT="(A,I12,A,I12)") &
711 460 : "Total number of neighbor interactions for process", para_env%mepos, ":", &
712 920 : nneighbor
713 460 : CALL compress(string)
714 460 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,T2,A)") TRIM(string)
715 :
716 460 : 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 19326 : 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 19326 : INTEGER, ALLOCATABLE, DIMENSION(:) :: indj
738 19326 : INTEGER, DIMENSION(:), POINTER :: work
739 19326 : INTEGER, DIMENSION(:, :), POINTER :: list_copy
740 : TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
741 :
742 19326 : NULLIFY (neighbor_kind_pair)
743 19326 : CALL timeset(routineN, handle)
744 : ! define a lookup table to get jkind for a given id_kind
745 57978 : ALLOCATE (indj(nkinds*(nkinds + 1)/2))
746 19326 : id_kind = 0
747 79605 : DO jkind = 1, nkinds
748 1004169 : DO ikind = jkind, nkinds
749 924564 : id_kind = id_kind + 1
750 984843 : indj(id_kind) = jkind
751 : END DO
752 : END DO
753 : ! loop over all nlists and sort the pairs within each list.
754 877886 : DO iab = 1, nonbonded%nlists
755 858560 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(iab)
756 858560 : npairs = neighbor_kind_pair%npairs
757 858560 : nscale = neighbor_kind_pair%nscale
758 858560 : IF (npairs /= 0) THEN
759 193417 : 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 574512 : ALLOCATE (work(1:npairs - nscale))
765 574512 : ALLOCATE (list_copy(2, 1:npairs - nscale))
766 : ! Copy of the pair list is required to perform the permutation below
767 : ! correctly.
768 906571516 : list_copy = neighbor_kind_pair%list(:, nscale + 1:npairs)
769 191504 : 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 151222922 : DO ipair = nscale + 1, npairs
773 151031418 : tmp = work(ipair - nscale)
774 151031418 : neighbor_kind_pair%list(1, ipair) = list_copy(1, tmp)
775 151222922 : neighbor_kind_pair%list(2, ipair) = list_copy(2, tmp)
776 : END DO
777 191504 : DEALLOCATE (work)
778 191504 : 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 193417 : max_alloc_size = nkinds*(nkinds + 1)/2 + nscale
789 193417 : IF (ASSOCIATED(neighbor_kind_pair%grp_kind_start)) THEN
790 59769 : DEALLOCATE (neighbor_kind_pair%grp_kind_start)
791 : END IF
792 580251 : ALLOCATE (neighbor_kind_pair%grp_kind_start(max_alloc_size))
793 193417 : IF (ASSOCIATED(neighbor_kind_pair%grp_kind_end)) THEN
794 59769 : DEALLOCATE (neighbor_kind_pair%grp_kind_end)
795 : END IF
796 386834 : ALLOCATE (neighbor_kind_pair%grp_kind_end(max_alloc_size))
797 193417 : IF (ASSOCIATED(neighbor_kind_pair%ij_kind)) THEN
798 59769 : DEALLOCATE (neighbor_kind_pair%ij_kind)
799 : END IF
800 580251 : ALLOCATE (neighbor_kind_pair%ij_kind(2, max_alloc_size))
801 : ! Start the first interval.
802 193417 : ipair = 1
803 193417 : neighbor_kind_pair%ngrp_kind = 1
804 193417 : neighbor_kind_pair%grp_kind_start(neighbor_kind_pair%ngrp_kind) = ipair
805 : ! Get ikind and jkind corresponding to id_kind.
806 193417 : id_kind = neighbor_kind_pair%id_kind(ipair)
807 193417 : jkind = indj(id_kind)
808 193417 : tmp = nkinds - jkind
809 193417 : ikind = nkinds + id_kind - nkinds*(nkinds + 1)/2 + (tmp*(tmp + 1)/2)
810 193417 : neighbor_kind_pair%ij_kind(1, neighbor_kind_pair%ngrp_kind) = ikind
811 193417 : neighbor_kind_pair%ij_kind(2, neighbor_kind_pair%ngrp_kind) = jkind
812 : ! Define the remaining intervals.
813 151821095 : DO ipair = 2, npairs
814 151821095 : IF (neighbor_kind_pair%id_kind(ipair) /= neighbor_kind_pair%id_kind(ipair - 1)) THEN
815 1116612 : neighbor_kind_pair%grp_kind_end(neighbor_kind_pair%ngrp_kind) = ipair - 1
816 1116612 : neighbor_kind_pair%ngrp_kind = neighbor_kind_pair%ngrp_kind + 1
817 1116612 : neighbor_kind_pair%grp_kind_start(neighbor_kind_pair%ngrp_kind) = ipair
818 : ! Get ikind and jkind corresponding to id_kind.
819 1116612 : id_kind = neighbor_kind_pair%id_kind(ipair)
820 1116612 : jkind = indj(id_kind)
821 1116612 : tmp = nkinds - jkind
822 1116612 : ikind = nkinds + id_kind - nkinds*(nkinds + 1)/2 + (tmp*(tmp + 1)/2)
823 1116612 : neighbor_kind_pair%ij_kind(1, neighbor_kind_pair%ngrp_kind) = ikind
824 1116612 : neighbor_kind_pair%ij_kind(2, neighbor_kind_pair%ngrp_kind) = jkind
825 : END IF
826 : END DO
827 : ! Finish the last interval.
828 193417 : 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 193417 : CALL reallocate(neighbor_kind_pair%grp_kind_start, 1, neighbor_kind_pair%ngrp_kind)
832 193417 : CALL reallocate(neighbor_kind_pair%grp_kind_end, 1, neighbor_kind_pair%ngrp_kind)
833 193417 : CALL reallocate(neighbor_kind_pair%ij_kind, 1, 2, 1, neighbor_kind_pair%ngrp_kind)
834 : END IF
835 : ! Clean the memory..
836 877886 : DEALLOCATE (neighbor_kind_pair%id_kind)
837 : END DO
838 19326 : DEALLOCATE (indj)
839 19326 : CALL timestop(handle)
840 19326 : END SUBROUTINE sort_neighbor_lists
841 :
842 0 : END MODULE fist_neighbor_lists
|
