Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Linked-cell neighbour finder with a Verlet skin for the NNP descriptor.
10 : !> Owns the per-nnp cell-list cache (positions, image pool, head/next
11 : !> chain, bin geometry), rebuilds or rebins it as atoms move, and walks
12 : !> it to fill the per-element neighbour buffers. Each &NNP force_eval
13 : !> owns its cache, so the lifetime tracks nnp_env_release.
14 : !> \author Dhruv Sharma (ds2173@cam.ac.uk)
15 : !> \author Christoph Schran (christoph.schran@rub.de)
16 : !> \date 2026-05-21
17 : ! **************************************************************************************************
18 : MODULE nnp_cell_list
19 :
20 : USE cell_types, ONLY: cell_type,&
21 : pbc,&
22 : real_to_scaled
23 : USE kinds, ONLY: dp
24 : USE nnp_environment_types, ONLY: nnp_cell_list_cache_type,&
25 : nnp_env_get,&
26 : nnp_neighbor_type,&
27 : nnp_type
28 : USE nnp_neighbor_interface, ONLY: nnp_neigh_grp_grow
29 : #include "./base/base_uses.f90"
30 :
31 : IMPLICIT NONE
32 :
33 : PRIVATE
34 :
35 : ! Verlet skin: defaults to MIN(0.5 bohr, 0.1*cutoff), or the VERLET_SKIN
36 : ! keyword value (nnp%verlet_skin). The chain rebuilds once an atom drifts
37 : ! more than skin/2 from its rebuild-time position.
38 : REAL(KIND=dp), PRIVATE, PARAMETER :: default_verlet_skin = 0.5_dp
39 : REAL(KIND=dp), PRIVATE, PARAMETER :: rebuild_eps = 1.0E-12_dp
40 :
41 : PUBLIC :: nnp_prepare_cell_list_cache, &
42 : nnp_compute_neighbors_cell_list
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief Prepare the cell-list cache for the current force evaluation: wrap
48 : !> positions into the primary cell, then either rebuild the bin geometry
49 : !> (cell, atom-count, cutoff or large-drift change) or rebin the chain in
50 : !> place while the Verlet skin still holds. Call once per force evaluation,
51 : !> before any nnp_compute_neighbors_cell_list call.
52 : !> \param nnp NNP environment whose cell_list_cache is prepared from current positions.
53 : !> \author Dhruv Sharma (ds2173@cam.ac.uk)
54 : ! **************************************************************************************************
55 55492 : SUBROUTINE nnp_prepare_cell_list_cache(nnp)
56 :
57 : TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
58 :
59 : INTEGER :: handle, i
60 : LOGICAL :: rebuild
61 : REAL(KIND=dp) :: dr2_max, exact_cutoff, list_cutoff, skin
62 : REAL(KIND=dp), DIMENSION(3) :: dr
63 : TYPE(cell_type), POINTER :: cell
64 :
65 55492 : CALL timeset('nnp_acsf_cell_list_prepare', handle)
66 :
67 55492 : NULLIFY (cell)
68 55492 : CALL nnp_env_get(nnp_env=nnp, cell=cell)
69 :
70 55492 : exact_cutoff = nnp%max_cut
71 55492 : IF (exact_cutoff > 0.0_dp) THEN
72 55492 : IF (nnp%verlet_skin >= 0.0_dp) THEN
73 : ! VERLET_SKIN keyword value (bohr).
74 0 : skin = nnp%verlet_skin
75 : ELSE
76 : ! Auto heuristic: cap at the default, scale down for small cutoffs.
77 55492 : skin = MIN(default_verlet_skin, 0.1_dp*exact_cutoff)
78 : END IF
79 : ELSE
80 0 : skin = 0.0_dp
81 : END IF
82 55492 : list_cutoff = exact_cutoff + skin
83 :
84 55492 : CALL nnp_cell_list_ensure_coord_buffers(nnp%cell_list_cache, nnp%num_atoms)
85 :
86 : ASSOCIATE (cache => nnp%cell_list_cache)
87 340374 : DO i = 1, nnp%num_atoms
88 1139528 : cache%coord_primary(:, i) = pbc(nnp%coord(:, i), cell, .TRUE.)
89 : CALL real_to_scaled(cache%coord_scaled(:, i), &
90 340374 : cache%coord_primary(:, i), cell)
91 : END DO
92 :
93 55492 : rebuild = .NOT. cache%initialized
94 55492 : IF (.NOT. rebuild) rebuild = (cache%num_atoms /= nnp%num_atoms)
95 55492 : IF (.NOT. rebuild) rebuild = (ABS(cache%exact_cutoff - exact_cutoff) > rebuild_eps)
96 55475 : IF (.NOT. rebuild) rebuild = (ABS(cache%list_cutoff - list_cutoff) > rebuild_eps)
97 55492 : IF (.NOT. rebuild) rebuild = (ABS(cache%verlet_skin - skin) > rebuild_eps)
98 221900 : IF (.NOT. rebuild) rebuild = ANY(cache%perd /= cell%perd)
99 55475 : IF (.NOT. rebuild) rebuild = (cache%orthorhombic .NEQV. cell%orthorhombic)
100 774980 : IF (.NOT. rebuild) rebuild = ANY(ABS(cache%hmat - cell%hmat) > rebuild_eps)
101 774126 : IF (.NOT. rebuild) rebuild = ANY(ABS(cache%h_inv - cell%h_inv) > rebuild_eps)
102 :
103 55281 : IF ((.NOT. rebuild) .AND. (skin > 0.0_dp)) THEN
104 55281 : dr2_max = 0.0_dp
105 299065 : DO i = 1, nnp%num_atoms
106 975136 : dr(:) = cache%coord_primary(:, i) - cache%ref_coord_primary(:, i)
107 975136 : dr(:) = pbc(dr, cell)
108 1030417 : dr2_max = MAX(dr2_max, DOT_PRODUCT(dr, dr))
109 : END DO
110 55281 : IF (dr2_max > 0.25_dp*skin**2) rebuild = .TRUE.
111 : END IF
112 :
113 55492 : IF (rebuild) THEN
114 54318 : CALL nnp_build_cell_list_cache(cache, cell, exact_cutoff, list_cutoff, skin)
115 1174 : ELSE IF (cache%initialized) THEN
116 : ! On reuse the bin geometry is frozen, but the head/next chain
117 : ! is rebinned from current positions so the walk still finds every
118 : ! in-cutoff pair.
119 1174 : CALL nnp_cell_list_rebin_chain(cache)
120 : END IF
121 : END ASSOCIATE
122 :
123 55492 : CALL timestop(handle)
124 :
125 55492 : END SUBROUTINE nnp_prepare_cell_list_cache
126 :
127 : ! **************************************************************************************************
128 : !> \brief Ensure the persistent coord buffers match num_atoms.
129 : !> Only (re-)allocates on first call or when num_atoms changes.
130 : !> \param cache cell-list cache whose coordinate buffers will be (re-)sized
131 : !> \param num_atoms number of primary-cell atoms the buffers must hold
132 : ! **************************************************************************************************
133 55492 : SUBROUTINE nnp_cell_list_ensure_coord_buffers(cache, num_atoms)
134 :
135 : TYPE(nnp_cell_list_cache_type), INTENT(INOUT) :: cache
136 : INTEGER, INTENT(IN) :: num_atoms
137 :
138 55492 : IF (ALLOCATED(cache%coord_primary)) THEN
139 55475 : IF (SIZE(cache%coord_primary, 2) == num_atoms) RETURN
140 0 : DEALLOCATE (cache%coord_primary)
141 : END IF
142 17 : IF (ALLOCATED(cache%coord_scaled)) DEALLOCATE (cache%coord_scaled)
143 17 : IF (ALLOCATED(cache%ref_coord_primary)) DEALLOCATE (cache%ref_coord_primary)
144 :
145 51 : ALLOCATE (cache%coord_primary(3, num_atoms))
146 34 : ALLOCATE (cache%coord_scaled(3, num_atoms))
147 34 : ALLOCATE (cache%ref_coord_primary(3, num_atoms))
148 :
149 : ! Resizing implies stale bin geometry; force the next prepare to rebuild.
150 17 : cache%initialized = .FALSE.
151 17 : cache%num_atoms = num_atoms
152 :
153 : END SUBROUTINE nnp_cell_list_ensure_coord_buffers
154 :
155 : ! **************************************************************************************************
156 : !> \brief Build the Cartesian linked-cell cache from the current positions.
157 : !>
158 : !> Positions must already live in cache%coord_primary / coord_scaled --
159 : !> this routine reads them in place and resets ref_coord_primary for
160 : !> the Verlet skin check.
161 : !> \param cache cell-list cache to populate; coord_primary / coord_scaled must already be filled
162 : !> \param cell simulation cell providing hmat, h_inv and per-axis periodicity
163 : !> \param exact_cutoff largest per-pair SF cutoff used to size the exact image ring
164 : !> \param list_cutoff cell-list bin width (exact_cutoff + Verlet skin)
165 : !> \param skin Verlet skin allowance used by the bin-reuse displacement check
166 : ! **************************************************************************************************
167 54318 : SUBROUTINE nnp_build_cell_list_cache(cache, cell, exact_cutoff, list_cutoff, skin)
168 :
169 : TYPE(nnp_cell_list_cache_type), INTENT(INOUT) :: cache
170 : TYPE(cell_type), INTENT(IN), POINTER :: cell
171 : REAL(KIND=dp), INTENT(IN) :: exact_cutoff, list_cutoff, skin
172 :
173 : INTEGER :: i, img, n_images, nx, ny, nz, tx, ty, tz
174 : INTEGER, DIMENSION(3) :: bin_index, shift
175 : REAL(KIND=dp) :: padding, safe_cutoff
176 : REAL(KIND=dp), DIMENSION(3) :: max_ref, min_ref, range_xyz, ref_pos
177 :
178 54318 : CALL nnp_cell_list_release_cell_structure(cache)
179 :
180 54318 : cache%exact_cutoff = exact_cutoff
181 54318 : cache%list_cutoff = list_cutoff
182 54318 : cache%verlet_skin = skin
183 217272 : cache%perd = cell%perd
184 54318 : cache%orthorhombic = cell%orthorhombic
185 706134 : cache%hmat = cell%hmat
186 706134 : cache%h_inv = cell%h_inv
187 :
188 1101062 : cache%ref_coord_primary(:, :) = cache%coord_primary(:, :)
189 :
190 54318 : CALL nnp_compute_pbc_copies(cache%exact_pbc_copies, cell, exact_cutoff)
191 54318 : CALL nnp_compute_pbc_copies(cache%list_pbc_copies, cell, list_cutoff)
192 : ! The bin stencil clips out-of-range bins rather than wrapping, so periodic
193 : ! neighbours come only from pre-replicated images. Force at least one ring
194 : ! on each periodic axis (MAX, not +) so boxes with L > 2*cutoff keep their
195 : ! wrap-around images.
196 217272 : cache%image_copies = MAX(cache%list_pbc_copies, cell%perd)
197 : ! The exact image ring must fit inside the replicated pool; otherwise the
198 : ! matcher would accept a ghost that was never created and drop a real pair.
199 217272 : CPASSERT(ALL(cache%exact_pbc_copies <= cache%image_copies))
200 :
201 54318 : nx = 2*cache%image_copies(1) + 1
202 54318 : ny = 2*cache%image_copies(2) + 1
203 54318 : nz = 2*cache%image_copies(3) + 1
204 54318 : n_images = cache%num_atoms*nx*ny*nz
205 54318 : cache%n_images = n_images
206 :
207 162954 : ALLOCATE (cache%image_atom(n_images))
208 162954 : ALLOCATE (cache%image_shift(3, n_images))
209 162954 : ALLOCATE (cache%image_translation(3, n_images))
210 :
211 217272 : min_ref(:) = HUGE(1.0_dp)
212 217272 : max_ref(:) = -HUGE(1.0_dp)
213 : img = 0
214 316004 : DO i = 1, cache%num_atoms
215 668314 : DO tx = -cache%image_copies(1), cache%image_copies(1)
216 1238178 : DO ty = -cache%image_copies(2), cache%image_copies(2)
217 2416290 : DO tz = -cache%image_copies(3), cache%image_copies(3)
218 1439798 : img = img + 1
219 5759192 : shift = [tx, ty, tz]
220 1439798 : cache%image_atom(img) = i
221 5759192 : cache%image_shift(:, img) = shift(:)
222 27356162 : cache%image_translation(:, img) = MATMUL(cell%hmat, REAL(shift, KIND=dp))
223 5759192 : ref_pos(:) = cache%coord_primary(:, i) + cache%image_translation(:, img)
224 5759192 : min_ref(:) = MIN(min_ref(:), ref_pos(:))
225 6383374 : max_ref(:) = MAX(max_ref(:), ref_pos(:))
226 : END DO
227 : END DO
228 : END DO
229 : END DO
230 :
231 54318 : padding = 0.5_dp*skin + rebuild_eps
232 217272 : cache%lower(:) = min_ref(:) - padding
233 217272 : cache%upper(:) = max_ref(:) + padding
234 217272 : range_xyz(:) = cache%upper(:) - cache%lower(:)
235 :
236 : ! Coarse grid: bin width = list_cutoff (bin_span 1, 27-bin walk). A finer
237 : ! grid changes the order in which neighbours are appended, which reorders
238 : ! the descriptor summation, so it is not used.
239 54318 : safe_cutoff = MAX(list_cutoff, rebuild_eps)
240 217272 : DO i = 1, 3
241 217272 : IF (range_xyz(i) > rebuild_eps) THEN
242 162954 : cache%nbin(i) = MAX(1, CEILING(range_xyz(i)/safe_cutoff))
243 162954 : cache%bin_width(i) = range_xyz(i)/REAL(cache%nbin(i), KIND=dp)
244 : cache%bin_span(i) = MIN(cache%nbin(i) - 1, &
245 162954 : CEILING(list_cutoff/MAX(cache%bin_width(i), rebuild_eps)))
246 : ELSE
247 0 : cache%nbin(i) = 1
248 0 : cache%bin_width(i) = 1.0_dp
249 0 : cache%bin_span(i) = 0
250 : END IF
251 : END DO
252 :
253 217272 : cache%n_cells = PRODUCT(cache%nbin)
254 162954 : ALLOCATE (cache%head(cache%n_cells))
255 108636 : ALLOCATE (cache%next(n_images))
256 215334 : cache%head(:) = 0
257 1494116 : cache%next(:) = 0
258 :
259 1494116 : DO img = 1, n_images
260 1439798 : i = cache%image_atom(img)
261 5759192 : ref_pos(:) = cache%coord_primary(:, i) + cache%image_translation(:, img)
262 1439798 : CALL nnp_cell_list_bin_from_position(cache, ref_pos, bin_index)
263 1439798 : cache%next(img) = cache%head(nnp_cell_list_linear_index(cache, bin_index))
264 1494116 : cache%head(nnp_cell_list_linear_index(cache, bin_index)) = img
265 : END DO
266 :
267 54318 : cache%initialized = .TRUE.
268 :
269 54318 : END SUBROUTINE nnp_build_cell_list_cache
270 :
271 : ! **************************************************************************************************
272 : !> \brief Rebin the head/next chain from current coord_primary while keeping the
273 : !> image pool, bin geometry and ref_coord_primary frozen. Used on Verlet
274 : !> skin reuse so the chain reflects atoms that crossed bin boundaries.
275 : !> \param cache cell-list cache whose head/next chain is rebuilt from current coord_primary
276 : ! **************************************************************************************************
277 1174 : SUBROUTINE nnp_cell_list_rebin_chain(cache)
278 :
279 : TYPE(nnp_cell_list_cache_type), INTENT(INOUT) :: cache
280 :
281 : INTEGER :: i, img, lin_idx
282 : INTEGER, DIMENSION(3) :: bin_index
283 : REAL(KIND=dp), DIMENSION(3) :: ref_pos
284 :
285 24386 : cache%head(:) = 0
286 513586 : cache%next(:) = 0
287 :
288 513586 : DO img = 1, cache%n_images
289 512412 : i = cache%image_atom(img)
290 2049648 : ref_pos(:) = cache%coord_primary(:, i) + cache%image_translation(:, img)
291 512412 : CALL nnp_cell_list_bin_from_position(cache, ref_pos, bin_index)
292 512412 : lin_idx = nnp_cell_list_linear_index(cache, bin_index)
293 512412 : cache%next(img) = cache%head(lin_idx)
294 513586 : cache%head(lin_idx) = img
295 : END DO
296 :
297 1174 : END SUBROUTINE nnp_cell_list_rebin_chain
298 :
299 : ! **************************************************************************************************
300 : !> \brief Free only the cell-structure pieces of the cache (bins and image
301 : !> metadata). The coordinate buffers are persistent and kept, so a
302 : !> rebuild can still read the current positions.
303 : !> \param cache cell-list cache whose bin geometry and image metadata are freed (coord buffers kept)
304 : ! **************************************************************************************************
305 54318 : SUBROUTINE nnp_cell_list_release_cell_structure(cache)
306 :
307 : TYPE(nnp_cell_list_cache_type), INTENT(INOUT) :: cache
308 :
309 54318 : IF (ALLOCATED(cache%image_atom)) DEALLOCATE (cache%image_atom)
310 54318 : IF (ALLOCATED(cache%head)) DEALLOCATE (cache%head)
311 54318 : IF (ALLOCATED(cache%next)) DEALLOCATE (cache%next)
312 54318 : IF (ALLOCATED(cache%image_shift)) DEALLOCATE (cache%image_shift)
313 54318 : IF (ALLOCATED(cache%image_translation)) DEALLOCATE (cache%image_translation)
314 :
315 54318 : cache%n_images = 0
316 54318 : cache%n_cells = 1
317 217272 : cache%exact_pbc_copies = 0
318 217272 : cache%list_pbc_copies = 0
319 217272 : cache%image_copies = 0
320 217272 : cache%nbin = 1
321 217272 : cache%bin_span = 0
322 217272 : cache%lower = 0.0_dp
323 217272 : cache%upper = 0.0_dp
324 217272 : cache%bin_width = 1.0_dp
325 :
326 54318 : END SUBROUTINE nnp_cell_list_release_cell_structure
327 :
328 : ! **************************************************************************************************
329 : !> \brief Fill the ACSF neighbour buffers for one central atom by walking the
330 : !> linked-cell stencil around its bin, pushing (j, dr, r) into the
331 : !> per-species-pair slabs for entries inside each pair cutoff. The own-bin
332 : !> self image (i == j, zero shift) is skipped.
333 : !> \param nnp NNP environment with cell_list_cache and neighbor_interface_state ready
334 : !> \param neighbor per-atom neighbour view filled in place (counters zeroed before this call)
335 : !> \param i central-atom index in the global ordering
336 : !> \author Dhruv Sharma (ds2173@cam.ac.uk)
337 : !> \author Christoph Schran (christoph.schran@rub.de)
338 : ! **************************************************************************************************
339 252785 : SUBROUTINE nnp_compute_neighbors_cell_list(nnp, neighbor, i)
340 :
341 : TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
342 : TYPE(nnp_neighbor_type), INTENT(INOUT) :: neighbor
343 : INTEGER, INTENT(IN) :: i
344 :
345 : INTEGER :: bx, by, bz, img, ind, j, neighbor_ind, &
346 : pair_slot, s
347 : INTEGER, DIMENSION(3) :: bin_index, current_shift
348 : REAL(KIND=dp) :: norm
349 : REAL(KIND=dp), DIMENSION(3) :: center, dr, image_pos
350 :
351 : ASSOCIATE (cache => nnp%cell_list_cache, &
352 : state => nnp%neighbor_interface_state)
353 :
354 252785 : ind = nnp%ele_ind(i)
355 1011140 : center(:) = cache%coord_primary(:, i)
356 252785 : CALL nnp_cell_list_bin_from_position(cache, center, bin_index)
357 :
358 684326 : DO bx = MAX(1, bin_index(1) - cache%bin_span(1)), &
359 505570 : MIN(cache%nbin(1), bin_index(1) + cache%bin_span(1))
360 1640194 : DO by = MAX(1, bin_index(2) - cache%bin_span(2)), &
361 684326 : MIN(cache%nbin(2), bin_index(2) + cache%bin_span(2))
362 5664609 : DO bz = MAX(1, bin_index(3) - cache%bin_span(3)), &
363 1640194 : MIN(cache%nbin(3), bin_index(3) + cache%bin_span(3))
364 17823824 : img = cache%head(nnp_cell_list_linear_index(cache, [bx, by, bz]))
365 103387781 : DO WHILE (img > 0)
366 97723172 : j = cache%image_atom(img)
367 390892688 : current_shift(:) = cache%image_shift(:, img)
368 171160725 : IF (j == i .AND. ALL(current_shift == 0)) THEN
369 252785 : img = cache%next(img)
370 252785 : CYCLE
371 : END IF
372 97470387 : IF (.NOT. nnp_cell_list_exact_image_match(cache, i, j, current_shift)) THEN
373 5102671 : img = cache%next(img)
374 5102671 : CYCLE
375 : END IF
376 :
377 369470864 : image_pos(:) = cache%coord_primary(:, j) + cache%image_translation(:, img)
378 369470864 : dr(:) = center(:) - image_pos(:)
379 369470864 : norm = NORM2(dr(:))
380 92367716 : neighbor_ind = nnp%ele_ind(j)
381 :
382 92367716 : IF (norm < state%pair_map(ind, neighbor_ind)%max_relevant_cutoff) THEN
383 7561404 : DO pair_slot = 1, state%pair_map(ind, neighbor_ind)%n_rad
384 3780702 : s = state%pair_map(ind, neighbor_ind)%rad_groups(pair_slot)
385 7561404 : IF (norm < nnp%rad(ind)%symfgrp(s)%cutoff) THEN
386 3780702 : neighbor%n_rad(s) = neighbor%n_rad(s) + 1
387 3780702 : IF (neighbor%n_rad(s) > neighbor%rad(s)%cap) THEN
388 199 : CALL nnp_neigh_grp_grow(neighbor%rad(s), neighbor%n_rad(s))
389 : END IF
390 3780702 : neighbor%rad(s)%ind(neighbor%n_rad(s)) = j
391 3780702 : neighbor%rad(s)%dist(1, neighbor%n_rad(s)) = dr(1)
392 3780702 : neighbor%rad(s)%dist(2, neighbor%n_rad(s)) = dr(2)
393 3780702 : neighbor%rad(s)%dist(3, neighbor%n_rad(s)) = dr(3)
394 3780702 : neighbor%rad(s)%dist(4, neighbor%n_rad(s)) = norm
395 : END IF
396 : END DO
397 :
398 8423152 : DO pair_slot = 1, state%pair_map(ind, neighbor_ind)%n_ang1
399 4642450 : s = state%pair_map(ind, neighbor_ind)%ang1_groups(pair_slot)
400 8423152 : IF (norm < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
401 4642450 : neighbor%n_ang1(s) = neighbor%n_ang1(s) + 1
402 4642450 : IF (neighbor%n_ang1(s) > neighbor%ang1(s)%cap) THEN
403 234 : CALL nnp_neigh_grp_grow(neighbor%ang1(s), neighbor%n_ang1(s))
404 : END IF
405 4642450 : neighbor%ang1(s)%ind(neighbor%n_ang1(s)) = j
406 4642450 : neighbor%ang1(s)%dist(1, neighbor%n_ang1(s)) = dr(1)
407 4642450 : neighbor%ang1(s)%dist(2, neighbor%n_ang1(s)) = dr(2)
408 4642450 : neighbor%ang1(s)%dist(3, neighbor%n_ang1(s)) = dr(3)
409 4642450 : neighbor%ang1(s)%dist(4, neighbor%n_ang1(s)) = norm
410 : END IF
411 : END DO
412 :
413 7071930 : DO pair_slot = 1, state%pair_map(ind, neighbor_ind)%n_ang2
414 3291228 : s = state%pair_map(ind, neighbor_ind)%ang2_groups(pair_slot)
415 7071930 : IF (norm < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
416 3291228 : neighbor%n_ang2(s) = neighbor%n_ang2(s) + 1
417 3291228 : IF (neighbor%n_ang2(s) > neighbor%ang2(s)%cap) THEN
418 212 : CALL nnp_neigh_grp_grow(neighbor%ang2(s), neighbor%n_ang2(s))
419 : END IF
420 3291228 : neighbor%ang2(s)%ind(neighbor%n_ang2(s)) = j
421 3291228 : neighbor%ang2(s)%dist(1, neighbor%n_ang2(s)) = dr(1)
422 3291228 : neighbor%ang2(s)%dist(2, neighbor%n_ang2(s)) = dr(2)
423 3291228 : neighbor%ang2(s)%dist(3, neighbor%n_ang2(s)) = dr(3)
424 3291228 : neighbor%ang2(s)%dist(4, neighbor%n_ang2(s)) = norm
425 : END IF
426 : END DO
427 : END IF
428 :
429 92367716 : img = cache%next(img)
430 : END DO
431 : END DO
432 : END DO
433 : END DO
434 :
435 : END ASSOCIATE
436 :
437 252785 : END SUBROUTINE nnp_compute_neighbors_cell_list
438 :
439 : ! **************************************************************************************************
440 : !> \brief Check whether a ghost image corresponds to one of the exact pbc copies.
441 : !> \param cache cell-list cache providing per-axis periodicity and scaled coordinates
442 : !> \param i central-atom index in the primary cell
443 : !> \param j candidate neighbour atom index in the primary cell
444 : !> \param image_shift 3-vector of integer cell translations identifying the ghost image
445 : !> \return .TRUE. if the (i, j, image_shift) image falls inside the exact_pbc_copies ring
446 : ! **************************************************************************************************
447 97470387 : LOGICAL FUNCTION nnp_cell_list_exact_image_match(cache, i, j, image_shift)
448 :
449 : TYPE(nnp_cell_list_cache_type), INTENT(IN) :: cache
450 : INTEGER, INTENT(IN) :: i, j
451 : INTEGER, DIMENSION(3), INTENT(IN) :: image_shift
452 :
453 : INTEGER :: d, minimal_shift
454 :
455 97470387 : nnp_cell_list_exact_image_match = .TRUE.
456 376798135 : DO d = 1, 3
457 376798135 : IF (cache%perd(d) == 1) THEN
458 282444029 : minimal_shift = NINT(cache%coord_scaled(d, i) - cache%coord_scaled(d, j))
459 282444029 : IF (ABS(minimal_shift - image_shift(d)) > cache%exact_pbc_copies(d)) THEN
460 97470387 : nnp_cell_list_exact_image_match = .FALSE.
461 : RETURN
462 : END IF
463 1986390 : ELSE IF (image_shift(d) /= 0) THEN
464 97470387 : nnp_cell_list_exact_image_match = .FALSE.
465 : RETURN
466 : END IF
467 : END DO
468 :
469 : END FUNCTION nnp_cell_list_exact_image_match
470 :
471 : ! **************************************************************************************************
472 : !> \brief Convert Cartesian position to linked-cell bin coordinates.
473 : !> \param cache cell-list cache providing the bin origin and bin widths
474 : !> \param pos Cartesian position (Bohr) to bin
475 : !> \param bin_index (out) 1-based bin indices clipped to [1, nbin] on each axis
476 : ! **************************************************************************************************
477 2204995 : SUBROUTINE nnp_cell_list_bin_from_position(cache, pos, bin_index)
478 :
479 : TYPE(nnp_cell_list_cache_type), INTENT(IN) :: cache
480 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: pos
481 : INTEGER, DIMENSION(3), INTENT(OUT) :: bin_index
482 :
483 : INTEGER :: d
484 : REAL(KIND=dp) :: scaled
485 :
486 8819980 : DO d = 1, 3
487 6614985 : scaled = (pos(d) - cache%lower(d))/MAX(cache%bin_width(d), rebuild_eps)
488 6614985 : bin_index(d) = INT(scaled) + 1
489 8819980 : bin_index(d) = MAX(1, MIN(cache%nbin(d), bin_index(d)))
490 : END DO
491 :
492 2204995 : END SUBROUTINE nnp_cell_list_bin_from_position
493 :
494 : ! **************************************************************************************************
495 : !> \brief Flatten 3D bin indices into the head/next storage index.
496 : !> \param cache cell-list cache providing the per-axis bin counts
497 : !> \param bin_index 1-based 3D bin indices to flatten
498 : !> \return linear index into cache%head
499 : ! **************************************************************************************************
500 6408166 : INTEGER FUNCTION nnp_cell_list_linear_index(cache, bin_index)
501 :
502 : TYPE(nnp_cell_list_cache_type), INTENT(IN) :: cache
503 : INTEGER, DIMENSION(3), INTENT(IN) :: bin_index
504 :
505 : nnp_cell_list_linear_index = bin_index(1) + cache%nbin(1)* &
506 6408166 : ((bin_index(2) - 1) + cache%nbin(2)*(bin_index(3) - 1))
507 :
508 6408166 : END FUNCTION nnp_cell_list_linear_index
509 :
510 : ! **************************************************************************************************
511 : !> \brief Number of PBC image rings along each lattice axis needed to cover a
512 : !> sphere of radius `cutoff` from any primary-cell atom. Uses the
513 : !> cell-vector projections onto the plane normals (correct for triclinic
514 : !> cells) and multiplies by cell%perd so non-periodic axes return 0.
515 : !> \param pbc_copies output: number of image rings on each of the 3 lattice axes
516 : !> \param cell cell whose hmat/h_inv/perd/deth are read
517 : !> \param cutoff sphere radius (Bohr) to cover
518 : !> \author Christoph Schran (christoph.schran@rub.de)
519 : ! **************************************************************************************************
520 108636 : SUBROUTINE nnp_compute_pbc_copies(pbc_copies, cell, cutoff)
521 :
522 : INTEGER, DIMENSION(3), INTENT(OUT) :: pbc_copies
523 : TYPE(cell_type), INTENT(IN), POINTER :: cell
524 : REAL(KIND=dp), INTENT(IN) :: cutoff
525 :
526 : REAL(KIND=dp) :: proja, projb, projc
527 : REAL(KIND=dp), DIMENSION(3) :: axb, axc, bxc
528 :
529 : ! A degenerate cell (deth <= 0) makes the NORM2 normalisations below
530 : ! divide by zero and the DO WHILE projection loops never terminate.
531 108636 : CPASSERT(cell%deth > 0.0_dp)
532 :
533 108636 : axb(1) = cell%hmat(2, 1)*cell%hmat(3, 2) - cell%hmat(3, 1)*cell%hmat(2, 2)
534 108636 : axb(2) = cell%hmat(3, 1)*cell%hmat(1, 2) - cell%hmat(1, 1)*cell%hmat(3, 2)
535 108636 : axb(3) = cell%hmat(1, 1)*cell%hmat(2, 2) - cell%hmat(2, 1)*cell%hmat(1, 2)
536 760452 : axb(:) = axb(:)/NORM2(axb(:))
537 :
538 108636 : axc(1) = cell%hmat(2, 1)*cell%hmat(3, 3) - cell%hmat(3, 1)*cell%hmat(2, 3)
539 108636 : axc(2) = cell%hmat(3, 1)*cell%hmat(1, 3) - cell%hmat(1, 1)*cell%hmat(3, 3)
540 108636 : axc(3) = cell%hmat(1, 1)*cell%hmat(2, 3) - cell%hmat(2, 1)*cell%hmat(1, 3)
541 760452 : axc(:) = axc(:)/NORM2(axc(:))
542 :
543 108636 : bxc(1) = cell%hmat(2, 2)*cell%hmat(3, 3) - cell%hmat(3, 2)*cell%hmat(2, 3)
544 108636 : bxc(2) = cell%hmat(3, 2)*cell%hmat(1, 3) - cell%hmat(1, 2)*cell%hmat(3, 3)
545 108636 : bxc(3) = cell%hmat(1, 2)*cell%hmat(2, 3) - cell%hmat(2, 2)*cell%hmat(1, 3)
546 760452 : bxc(:) = bxc(:)/NORM2(bxc(:))
547 :
548 434544 : proja = ABS(SUM(cell%hmat(:, 1)*bxc(:)))*0.5_dp
549 434544 : projb = ABS(SUM(cell%hmat(:, 2)*axc(:)))*0.5_dp
550 434544 : projc = ABS(SUM(cell%hmat(:, 3)*axb(:)))*0.5_dp
551 :
552 108636 : pbc_copies(:) = 0
553 433632 : DO WHILE ((pbc_copies(1) + 1)*proja <= cutoff)
554 324996 : pbc_copies(1) = pbc_copies(1) + 1
555 : END DO
556 433632 : DO WHILE ((pbc_copies(2) + 1)*projb <= cutoff)
557 324996 : pbc_copies(2) = pbc_copies(2) + 1
558 : END DO
559 433640 : DO WHILE ((pbc_copies(3) + 1)*projc <= cutoff)
560 325004 : pbc_copies(3) = pbc_copies(3) + 1
561 : END DO
562 434544 : pbc_copies(:) = pbc_copies(:)*cell%perd(:)
563 :
564 108636 : END SUBROUTINE nnp_compute_pbc_copies
565 :
566 : END MODULE nnp_cell_list
|
