(git:98357aa)
Loading...
Searching...
No Matches
nnp_acsf.F
Go to the documentation of this file.
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 Functionality for atom centered symmetry functions
10!> for neural network potentials
11!> \author Christoph Schran (christoph.schran@rub.de)
12!> \author Dhruv Sharma (ds2173@cam.ac.uk)
13!> \date 2020-10-10
14! **************************************************************************************************
15MODULE nnp_acsf
16 USE cp_log_handling, ONLY: cp_get_default_logger,&
17 cp_logger_get_default_unit_nr,&
18 cp_logger_type
19 USE kinds, ONLY: default_string_length,&
20 dp
21 USE mathconstants, ONLY: pi
22 USE message_passing, ONLY: mp_para_env_type
23 USE nnp_cell_list, ONLY: nnp_compute_neighbors_cell_list,&
24 nnp_prepare_cell_list_cache
25 USE nnp_environment_types, ONLY: nnp_acsf_ang_type,&
26 nnp_acsf_rad_type,&
27 nnp_cut_cos,&
28 nnp_cut_tanh,&
29 nnp_symfgrp_type,&
30 nnp_type
31 USE nnp_neighbor_interface, ONLY: nnp_grp_grow_dgdr,&
32 nnp_neighbor_interface_prepare,&
33 nnp_neighbor_interface_reset_neighbor,&
34 nnp_workspace_grow_caches
35 USE periodic_table, ONLY: get_ptable_info
36
37!$ USE omp_lib, ONLY: omp_get_max_threads
38#include "./base/base_uses.f90"
39
40 IMPLICIT NONE
41
42 PRIVATE
43
44 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .false.
45 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nnp_acsf'
46
47 ! Cutoff-equality tolerance for grouping symmetry functions in
48 ! nnp_init_acsf_groups: SFs whose cutoffs agree to within this absolute
49 ! tolerance share a group (and a spline grid).
50 REAL(KIND=dp), PARAMETER, PRIVATE :: cutoff_eq_tol = 1.0e-5_dp
51
52 ! Public subroutines ***
53 PUBLIC :: nnp_calc_acsf, &
54 nnp_prepare_neighbor_cache, &
55 nnp_init_acsf_groups, &
56 nnp_sort_acsf, &
57 nnp_sort_ele, &
58 nnp_write_acsf
59
60CONTAINS
61
62! **************************************************************************************************
63!> \brief Calculate atom centered symmetry functions for given atom i
64!>
65!> Per-atom symmetry-function gradients live in the per-element neighbour
66!> workspace as sparse per-group arrays: self_dGdr (atom i), dGdr_rad (radial
67!> group s), and dGdr_ang_jj / dGdr_ang_kk (angular group s, j- and k-side),
68!> each addressed by the neighbour's slot in workspace(ind)%neighbor. No global
69!> (3, n_sf, num_atoms) slab is needed: only atoms in atom i's neighbour lists
70!> get derivatives.
71!>
72!> \param nnp NNP environment with persistent neighbour caches populated by
73!> nnp_prepare_neighbor_cache (must be called once before the per-atom loop).
74!> \param i central-atom index in the global atom ordering; selects
75!> nnp%ele_ind(i) for per-element scratch routing.
76!> \param calc_forces if .TRUE., populate the per-element dGdr workspaces for caller-side
77!> force assembly via nnp_scatter_dgdr_to_forces.
78!> \param stress optional per-input-node stress accumulator (only valid when calc_forces).
79!> \date 2020-10-10
80!> \author Christoph Schran (christoph.schran@rub.de)
81!> \author Dhruv Sharma (ds2173@cam.ac.uk)
82! **************************************************************************************************
83 SUBROUTINE nnp_calc_acsf(nnp, i, calc_forces, stress)
84
85 TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
86 INTEGER, INTENT(IN) :: i
87 LOGICAL, INTENT(IN) :: calc_forces
88 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
89 OPTIONAL :: stress
90
91 CHARACTER(len=*), PARAMETER :: routinen = 'nnp_calc_acsf'
92
93 INTEGER :: handle, handle_nlist, handle_sf, ii, ind, izeta_il, j, k, l, m, n_ang1_s, &
94 n_ang2_s, n_input_nodes, n_symf_s, nthreads_ang, off, peak, s, sf
95 LOGICAL :: do_forces, homo_grp
96 REAL(kind=dp) :: angular_il, arg_il, costheta_il, cutoff_s, cutoff_sqr, dfcut3_il, &
97 dfcutdr1_il, dfcutdr2_il, dfcutdr3_il, dgdx_t1, dgdx_t2, dsymdr1_il, dsymdr2_il, &
98 dsymdr3_il, eta_il, f_il, fcut3_il, ftot_il, g_il, inv_g2_il, lam_il, pref_il, &
99 pref_lam_il, prefzeta_il, r1, r1_inv, r2, r2_inv, r2sum_il, r3, r3_inv, r3_sqr, rsqr1, &
100 rsqr2, rsqr3, sym_il, symtmp_il, tanh_il, tmp1_il, tmp2_il, tmp3_il, tmp_il, tmpzeta_il, &
101 zeta_il
102 REAL(kind=dp), DIMENSION(3) :: dcosbase1_il, dcosbase2_il, &
103 dcosbase3_il, dr1dx_il, dr2dx_il, &
104 dr3dx_il, f_jj_il, f_kk_il, rvect1, &
105 rvect2, rvect3
106
107! Inlined angular kernel variables (eliminates nnp_calc_ang call overhead
108! and intermediate angular_force3tmp array on the serial force path).
109
110 CALL timeset(routinen, handle)
111
112 !determine index of atom type
113 ind = nnp%ele_ind(i)
114 do_forces = calc_forces
115
116 ! Lazy one-shot build of the per-radial-group spline tables. They depend
117 ! only on eta/rs/cutoff/cut_type, fixed once nnp_init_acsf_groups and
118 ! nnp_sort_acsf have run; %spline_built is per-nnp so re-init is safe.
119 ! Guard the indexed access in case this element has no radial groups.
120 IF (nnp%rad(ind)%n_symfgrp > 0) THEN
121 IF (.NOT. nnp%rad(ind)%symfgrp(1)%spline_built) CALL nnp_build_radial_splines(nnp)
122 END IF
123
124 ! Persistent per-element workspace bindings. The dGdr_* slabs grow lazily
125 ! via nnp_grp_grow_dGdr; the fc_cache*/dfc_cache* buffers are bound in an
126 ! inner ASSOCIATE after nnp_workspace_grow_caches has sized them, so the
127 ! alias stays valid.
128 associate(workspace => nnp%neighbor_interface_state%workspace(ind), &
129 neighbor => nnp%neighbor_interface_state%workspace(ind)%neighbor, &
130 radial_symtmp => nnp%neighbor_interface_state%workspace(ind)%radial_sym, &
131 radial_forcetmp => nnp%neighbor_interface_state%workspace(ind)%radial_force, &
132 angular_symtmp => nnp%neighbor_interface_state%workspace(ind)%angular_sym, &
133 angular_force3tmp => nnp%neighbor_interface_state%workspace(ind)%angular_force, &
134 self_dgdr => nnp%neighbor_interface_state%workspace(ind)%self_dGdr)
135
136 n_input_nodes = nnp%neighbor_interface_state%workspace(ind)%n_input_nodes
137 IF (do_forces) self_dgdr(:, 1:n_input_nodes) = 0.0_dp
138
139 ! Walk the linked-cell candidates directly. The cell-list/neighbour cache
140 ! is prepared once per force evaluation by nnp_prepare_neighbor_cache, not
141 ! here, so the per-atom walk stays O(neighbours).
142 CALL timeset('nnp_acsf_neighbor_fill', handle_nlist)
143 neighbor%pbc_copies = nnp%cell_list_cache%exact_pbc_copies
144 CALL nnp_neighbor_interface_reset_neighbor(nnp, ind)
145 CALL nnp_compute_neighbors_cell_list(nnp, neighbor, i)
146 CALL timestop(handle_nlist)
147
148 ! Reset y:
149 nnp%rad(ind)%y = 0.0_dp
150 nnp%ang(ind)%y = 0.0_dp
151
152 ! Grow the per-element 1D angular cutoff caches to this atom's peak
153 ! per-group neighbour count, then bind them in an inner ASSOCIATE.
154 ! MAXVAL of a zero-size array is compiler-defined, so guard each axis.
155 peak = 0
156 IF (SIZE(neighbor%n_ang1) > 0) peak = max(peak, maxval(neighbor%n_ang1))
157 IF (SIZE(neighbor%n_ang2) > 0) peak = max(peak, maxval(neighbor%n_ang2))
158 IF (peak > 0) CALL nnp_workspace_grow_caches(workspace, peak)
159
160 associate(fc_cache1 => workspace%fc_cache1, &
161 dfc_cache1 => workspace%dfc_cache1, &
162 fc_cache2 => workspace%fc_cache2, &
163 dfc_cache2 => workspace%dfc_cache2)
164
165 !calc forces
166 IF (do_forces) THEN
167 !loop over radial sym fnct grps
168 CALL timeset('nnp_acsf_radial', handle_sf)
169 DO s = 1, nnp%rad(ind)%n_symfgrp
170 n_symf_s = nnp%rad(ind)%symfgrp(s)%n_symf
171 ! Per-group dense buffer: (3, n_symf_s, cap_s). Grown lazily.
172 CALL nnp_grp_grow_dgdr(workspace%dGdr_rad(s), neighbor%n_rad(s))
173 associate(rad_buf => workspace%dGdr_rad(s)%data)
174 !loop over associated neighbors
175 DO j = 1, neighbor%n_rad(s)
176 rvect1 = neighbor%rad(s)%dist(1:3, j)
177 r1 = neighbor%rad(s)%dist(4, j)
178 CALL nnp_calc_rad(nnp, ind, s, rvect1, r1, &
179 radial_symtmp(1:n_symf_s), &
180 radial_forcetmp(:, 1:n_symf_s))
181 ! Per-group dense write: rad_buf(:, sf, j) holds dG_m/dr_j.
182 DO sf = 1, n_symf_s
183 m = nnp%rad(ind)%symfgrp(s)%symf(sf)
184 self_dgdr(1, m) = self_dgdr(1, m) + radial_forcetmp(1, sf)
185 self_dgdr(2, m) = self_dgdr(2, m) + radial_forcetmp(2, sf)
186 self_dgdr(3, m) = self_dgdr(3, m) + radial_forcetmp(3, sf)
187 rad_buf(1, sf, j) = -radial_forcetmp(1, sf)
188 rad_buf(2, sf, j) = -radial_forcetmp(2, sf)
189 rad_buf(3, sf, j) = -radial_forcetmp(3, sf)
190 IF (PRESENT(stress)) THEN
191 DO l = 1, 3
192 stress(:, l, m) = stress(:, l, m) + rvect1(:)*radial_forcetmp(l, sf)
193 END DO
194 END IF
195 nnp%rad(ind)%y(m) = nnp%rad(ind)%y(m) + radial_symtmp(sf)
196 END DO
197 END DO
198 END associate
199 END DO
200 CALL timestop(handle_sf)
201
202 !loop over angular sym fnct grps
203 CALL timeset('nnp_acsf_angular', handle_sf)
204 off = nnp%n_rad(ind)
205
206 ! OpenMP over the angular group index s, taken only with >1 thread
207 ! and >1 group. Angular groups partition the input-node index m
208 ! disjointly (nnp_sort_acsf), so self_dGdr, stress and y are written
209 ! without races. The serial path below is identical.
210 nthreads_ang = 1
211!$ nthreads_ang = omp_get_max_threads()
212 IF (nthreads_ang > 1 .AND. nnp%ang(ind)%n_symfgrp > 1) THEN
213 IF (PRESENT(stress)) THEN
214 CALL nnp_acsf_angular_loop_omp(nnp, ind, self_dgdr, off, stress)
215 ELSE
216 CALL nnp_acsf_angular_loop_omp(nnp, ind, self_dgdr, off)
217 END IF
218 ELSE
219 DO s = 1, nnp%ang(ind)%n_symfgrp
220 cutoff_s = nnp%ang(ind)%symfgrp(s)%cutoff
221 cutoff_sqr = cutoff_s*cutoff_s
222 n_symf_s = nnp%ang(ind)%symfgrp(s)%n_symf
223 n_ang1_s = neighbor%n_ang1(s)
224 homo_grp = (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2))
225
226 ! Grow per-group dense buffers. jj is always indexed in ang1.
227 ! kk is indexed in ang1 for homo groups and in ang2 for hetero.
228 CALL nnp_grp_grow_dgdr(workspace%dGdr_ang_jj(s), n_ang1_s)
229 IF (homo_grp) THEN
230 CALL nnp_grp_grow_dgdr(workspace%dGdr_ang_kk(s), n_ang1_s)
231 n_ang2_s = 0
232 ELSE
233 n_ang2_s = neighbor%n_ang2(s)
234 CALL nnp_grp_grow_dgdr(workspace%dGdr_ang_kk(s), n_ang2_s)
235 END IF
236
237 ! Per-group reset. Triplets accumulate into the same (sf, j) slot
238 ! across multiple k partners, so we MUST zero before the triplet loop.
239 IF (n_ang1_s > 0) workspace%dGdr_ang_jj(s)%data(:, 1:n_symf_s, 1:n_ang1_s) = 0.0_dp
240 IF (homo_grp) THEN
241 IF (n_ang1_s > 0) workspace%dGdr_ang_kk(s)%data(:, 1:n_symf_s, 1:n_ang1_s) = 0.0_dp
242 ELSE
243 IF (n_ang2_s > 0) workspace%dGdr_ang_kk(s)%data(:, 1:n_symf_s, 1:n_ang2_s) = 0.0_dp
244 END IF
245
246 ! Precompute cutoff values and derivatives for ang1 neighbors
247 CALL nnp_fill_fc_dfc_cache(neighbor%ang1(s)%dist, n_ang1_s, &
248 nnp%cut_type, cutoff_s, fc_cache1, dfc_cache1)
249
250 ! Inlined angular kernel: compute + scatter fused. The (j,k)
251 ! geometry is computed once, then the SF loop scatters sym
252 ! values and forces directly into self_dGdr / jj_buf / kk_buf,
253 ! keeping per-triple scalars in registers. Inlined equivalent
254 ! of nnp_calc_ang; the OMP path calls it directly.
255 IF (homo_grp) THEN
256 associate(jj_buf => workspace%dGdr_ang_jj(s)%data, &
257 kk_buf => workspace%dGdr_ang_kk(s)%data, &
258 grp_il => nnp%ang(ind)%symfgrp(s))
259 DO j = 1, n_ang1_s
260 rvect1 = neighbor%ang1(s)%dist(1:3, j)
261 r1 = neighbor%ang1(s)%dist(4, j)
262 DO k = j + 1, n_ang1_s
263 rvect2 = neighbor%ang1(s)%dist(1:3, k)
264 r2 = neighbor%ang1(s)%dist(4, k)
265 rvect3(1) = rvect2(1) - rvect1(1)
266 rvect3(2) = rvect2(2) - rvect1(2)
267 rvect3(3) = rvect2(3) - rvect1(3)
268 r3_sqr = rvect3(1)*rvect3(1) + rvect3(2)*rvect3(2) + rvect3(3)*rvect3(3)
269 IF (r3_sqr < cutoff_sqr) THEN
270 r3 = sqrt(r3_sqr)
271
272 ! -- per-triple geometry (once) --
273 rsqr1 = r1*r1; rsqr2 = r2*r2; rsqr3 = r3*r3
274 r2sum_il = rsqr1 + rsqr2 + rsqr3
275 f_il = rsqr3 - rsqr1 - rsqr2
276 g_il = -2.0_dp*r1*r2
277 costheta_il = f_il/g_il
278
279 SELECT CASE (nnp%cut_type)
280 CASE (nnp_cut_cos)
281 arg_il = pi*r3/cutoff_s
282 fcut3_il = 0.5_dp*(cos(arg_il) + 1.0_dp)
283 dfcut3_il = -0.5_dp*sin(arg_il)*(pi/cutoff_s)
284 CASE (nnp_cut_tanh)
285 tanh_il = tanh(1.0_dp - r3/cutoff_s)
286 fcut3_il = tanh_il**3
287 dfcut3_il = (-3.0_dp/cutoff_s)*(tanh_il**2 - tanh_il**4)
288 CASE DEFAULT
289 cpabort("NNP| Cutoff function unknown")
290 END SELECT
291
292 ftot_il = fc_cache1(j)*fc_cache1(k)*fcut3_il
293 dfcutdr1_il = dfc_cache1(j)*fc_cache1(k)*fcut3_il
294 dfcutdr2_il = fc_cache1(j)*dfc_cache1(k)*fcut3_il
295 dfcutdr3_il = fc_cache1(j)*fc_cache1(k)*dfcut3_il
296
297 r1_inv = 1.0_dp/r1; r2_inv = 1.0_dp/r2; r3_inv = 1.0_dp/r3
298 dr1dx_il(:) = rvect1(:)*r1_inv
299 dr2dx_il(:) = rvect2(:)*r2_inv
300 dr3dx_il(:) = rvect3(:)*r3_inv
301
302 inv_g2_il = 1.0_dp/(g_il*g_il)
303 DO ii = 1, 3
304 dgdx_t1 = 2.0_dp*r2*dr1dx_il(ii)
305 dgdx_t2 = 2.0_dp*r1*dr2dx_il(ii)
306 dcosbase1_il(ii) = -2.0_dp*(rvect1(ii) + rvect2(ii))*g_il &
307 - f_il*(-(dgdx_t1 + dgdx_t2))
308 dcosbase2_il(ii) = 2.0_dp*(rvect3(ii) + rvect1(ii))*g_il &
309 - f_il*dgdx_t1
310 dcosbase3_il(ii) = 2.0_dp*(rvect2(ii) - rvect3(ii))*g_il &
311 - f_il*dgdx_t2
312 END DO
313
314 ! -- fused SF loop: compute + direct scatter --
315 DO sf = 1, n_symf_s
316 m = off + grp_il%symf(sf)
317 lam_il = grp_il%pack_lam(sf)
318 zeta_il = grp_il%pack_zeta(sf)
319 eta_il = grp_il%pack_eta(sf)
320 prefzeta_il = grp_il%pack_prefzeta(sf)
321
322 tmp_il = 1.0_dp + lam_il*costheta_il
323 IF (tmp_il <= 0.0_dp) THEN
324 tmpzeta_il = 0.0_dp
325 angular_il = 0.0_dp
326 ELSE
327 IF (grp_il%pack_use_int_zeta(sf)) THEN
328 izeta_il = grp_il%pack_izeta(sf)
329 tmpzeta_il = tmp_il**(izeta_il - 1)
330 ELSE
331 tmpzeta_il = tmp_il**(zeta_il - 1.0_dp)
332 END IF
333 angular_il = tmpzeta_il*tmp_il
334 END IF
335
336 symtmp_il = exp(-eta_il*r2sum_il)
337 sym_il = prefzeta_il*angular_il*symtmp_il*ftot_il
338 nnp%ang(ind)%y(m - off) = nnp%ang(ind)%y(m - off) + sym_il
339
340 pref_lam_il = zeta_il*tmpzeta_il*lam_il*inv_g2_il
341 tmp_il = -2.0_dp*symtmp_il*eta_il
342 dsymdr1_il = tmp_il*r1
343 dsymdr2_il = tmp_il*r2
344 dsymdr3_il = tmp_il*r3
345
346 pref_il = prefzeta_il*symtmp_il*ftot_il
347 tmp1_il = prefzeta_il*angular_il*(ftot_il*dsymdr1_il + dfcutdr1_il*symtmp_il)
348 tmp2_il = prefzeta_il*angular_il*(ftot_il*dsymdr2_il + dfcutdr2_il*symtmp_il)
349 tmp3_il = prefzeta_il*angular_il*(ftot_il*dsymdr3_il + dfcutdr3_il*symtmp_il)
350
351 DO ii = 1, 3
352 f_jj_il(ii) = pref_il*pref_lam_il*dcosbase2_il(ii) &
353 - tmp1_il*dr1dx_il(ii) + tmp3_il*dr3dx_il(ii)
354 f_kk_il(ii) = pref_il*pref_lam_il*dcosbase3_il(ii) &
355 - tmp2_il*dr2dx_il(ii) - tmp3_il*dr3dx_il(ii)
356 self_dgdr(ii, m) = self_dgdr(ii, m) &
357 + pref_il*pref_lam_il*dcosbase1_il(ii) &
358 + tmp1_il*dr1dx_il(ii) + tmp2_il*dr2dx_il(ii)
359 jj_buf(ii, sf, j) = jj_buf(ii, sf, j) + f_jj_il(ii)
360 kk_buf(ii, sf, k) = kk_buf(ii, sf, k) + f_kk_il(ii)
361 END DO
362 IF (PRESENT(stress)) THEN
363 DO l = 1, 3
364 stress(:, l, m) = stress(:, l, m) &
365 - rvect1(:)*f_jj_il(l) - rvect2(:)*f_kk_il(l)
366 END DO
367 END IF
368 END DO
369
370 END IF
371 END DO
372 END DO
373 END associate
374 ELSE
375 ! Precompute cutoff values for ang2 neighbors (different elements)
376 CALL nnp_fill_fc_dfc_cache(neighbor%ang2(s)%dist, n_ang2_s, &
377 nnp%cut_type, cutoff_s, fc_cache2, dfc_cache2)
378
379 associate(jj_buf => workspace%dGdr_ang_jj(s)%data, &
380 kk_buf => workspace%dGdr_ang_kk(s)%data, &
381 grp_il => nnp%ang(ind)%symfgrp(s))
382 DO j = 1, n_ang1_s
383 rvect1 = neighbor%ang1(s)%dist(1:3, j)
384 r1 = neighbor%ang1(s)%dist(4, j)
385 DO k = 1, n_ang2_s
386 rvect2 = neighbor%ang2(s)%dist(1:3, k)
387 r2 = neighbor%ang2(s)%dist(4, k)
388 rvect3(1) = rvect2(1) - rvect1(1)
389 rvect3(2) = rvect2(2) - rvect1(2)
390 rvect3(3) = rvect2(3) - rvect1(3)
391 r3_sqr = rvect3(1)*rvect3(1) + rvect3(2)*rvect3(2) + rvect3(3)*rvect3(3)
392 IF (r3_sqr < cutoff_sqr) THEN
393 r3 = sqrt(r3_sqr)
394
395 ! -- per-triple geometry (once) --
396 rsqr1 = r1*r1; rsqr2 = r2*r2; rsqr3 = r3*r3
397 r2sum_il = rsqr1 + rsqr2 + rsqr3
398 f_il = rsqr3 - rsqr1 - rsqr2
399 g_il = -2.0_dp*r1*r2
400 costheta_il = f_il/g_il
401
402 SELECT CASE (nnp%cut_type)
403 CASE (nnp_cut_cos)
404 arg_il = pi*r3/cutoff_s
405 fcut3_il = 0.5_dp*(cos(arg_il) + 1.0_dp)
406 dfcut3_il = -0.5_dp*sin(arg_il)*(pi/cutoff_s)
407 CASE (nnp_cut_tanh)
408 tanh_il = tanh(1.0_dp - r3/cutoff_s)
409 fcut3_il = tanh_il**3
410 dfcut3_il = (-3.0_dp/cutoff_s)*(tanh_il**2 - tanh_il**4)
411 CASE DEFAULT
412 cpabort("NNP| Cutoff function unknown")
413 END SELECT
414
415 ftot_il = fc_cache1(j)*fc_cache2(k)*fcut3_il
416 dfcutdr1_il = dfc_cache1(j)*fc_cache2(k)*fcut3_il
417 dfcutdr2_il = fc_cache1(j)*dfc_cache2(k)*fcut3_il
418 dfcutdr3_il = fc_cache1(j)*fc_cache2(k)*dfcut3_il
419
420 r1_inv = 1.0_dp/r1; r2_inv = 1.0_dp/r2; r3_inv = 1.0_dp/r3
421 dr1dx_il(:) = rvect1(:)*r1_inv
422 dr2dx_il(:) = rvect2(:)*r2_inv
423 dr3dx_il(:) = rvect3(:)*r3_inv
424
425 inv_g2_il = 1.0_dp/(g_il*g_il)
426 DO ii = 1, 3
427 dgdx_t1 = 2.0_dp*r2*dr1dx_il(ii)
428 dgdx_t2 = 2.0_dp*r1*dr2dx_il(ii)
429 dcosbase1_il(ii) = -2.0_dp*(rvect1(ii) + rvect2(ii))*g_il &
430 - f_il*(-(dgdx_t1 + dgdx_t2))
431 dcosbase2_il(ii) = 2.0_dp*(rvect3(ii) + rvect1(ii))*g_il &
432 - f_il*dgdx_t1
433 dcosbase3_il(ii) = 2.0_dp*(rvect2(ii) - rvect3(ii))*g_il &
434 - f_il*dgdx_t2
435 END DO
436
437 ! -- fused SF loop: compute + direct scatter --
438 DO sf = 1, n_symf_s
439 m = off + grp_il%symf(sf)
440 lam_il = grp_il%pack_lam(sf)
441 zeta_il = grp_il%pack_zeta(sf)
442 eta_il = grp_il%pack_eta(sf)
443 prefzeta_il = grp_il%pack_prefzeta(sf)
444
445 tmp_il = 1.0_dp + lam_il*costheta_il
446 IF (tmp_il <= 0.0_dp) THEN
447 tmpzeta_il = 0.0_dp
448 angular_il = 0.0_dp
449 ELSE
450 IF (grp_il%pack_use_int_zeta(sf)) THEN
451 izeta_il = grp_il%pack_izeta(sf)
452 tmpzeta_il = tmp_il**(izeta_il - 1)
453 ELSE
454 tmpzeta_il = tmp_il**(zeta_il - 1.0_dp)
455 END IF
456 angular_il = tmpzeta_il*tmp_il
457 END IF
458
459 symtmp_il = exp(-eta_il*r2sum_il)
460 sym_il = prefzeta_il*angular_il*symtmp_il*ftot_il
461 nnp%ang(ind)%y(m - off) = nnp%ang(ind)%y(m - off) + sym_il
462
463 pref_lam_il = zeta_il*tmpzeta_il*lam_il*inv_g2_il
464 tmp_il = -2.0_dp*symtmp_il*eta_il
465 dsymdr1_il = tmp_il*r1
466 dsymdr2_il = tmp_il*r2
467 dsymdr3_il = tmp_il*r3
468
469 pref_il = prefzeta_il*symtmp_il*ftot_il
470 tmp1_il = prefzeta_il*angular_il*(ftot_il*dsymdr1_il + dfcutdr1_il*symtmp_il)
471 tmp2_il = prefzeta_il*angular_il*(ftot_il*dsymdr2_il + dfcutdr2_il*symtmp_il)
472 tmp3_il = prefzeta_il*angular_il*(ftot_il*dsymdr3_il + dfcutdr3_il*symtmp_il)
473
474 DO ii = 1, 3
475 f_jj_il(ii) = pref_il*pref_lam_il*dcosbase2_il(ii) &
476 - tmp1_il*dr1dx_il(ii) + tmp3_il*dr3dx_il(ii)
477 f_kk_il(ii) = pref_il*pref_lam_il*dcosbase3_il(ii) &
478 - tmp2_il*dr2dx_il(ii) - tmp3_il*dr3dx_il(ii)
479 self_dgdr(ii, m) = self_dgdr(ii, m) &
480 + pref_il*pref_lam_il*dcosbase1_il(ii) &
481 + tmp1_il*dr1dx_il(ii) + tmp2_il*dr2dx_il(ii)
482 jj_buf(ii, sf, j) = jj_buf(ii, sf, j) + f_jj_il(ii)
483 kk_buf(ii, sf, k) = kk_buf(ii, sf, k) + f_kk_il(ii)
484 END DO
485 IF (PRESENT(stress)) THEN
486 DO l = 1, 3
487 stress(:, l, m) = stress(:, l, m) &
488 - rvect1(:)*f_jj_il(l) - rvect2(:)*f_kk_il(l)
489 END DO
490 END IF
491 END DO
492
493 END IF
494 END DO
495 END DO
496 END associate
497 END IF
498 END DO
499 END IF
500 CALL timestop(handle_sf)
501 ELSE
502 !loop over radial sym fnct grps
503 CALL timeset('nnp_acsf_radial', handle_sf)
504 DO s = 1, nnp%rad(ind)%n_symfgrp
505 !loop over associated neighbors
506 DO j = 1, neighbor%n_rad(s)
507 rvect1 = neighbor%rad(s)%dist(1:3, j)
508 r1 = neighbor%rad(s)%dist(4, j)
509 CALL nnp_calc_rad(nnp, ind, s, rvect1, r1, radial_symtmp(1:nnp%rad(ind)%symfgrp(s)%n_symf))
510 DO sf = 1, nnp%rad(ind)%symfgrp(s)%n_symf
511 m = nnp%rad(ind)%symfgrp(s)%symf(sf)
512 nnp%rad(ind)%y(m) = nnp%rad(ind)%y(m) + radial_symtmp(sf)
513 END DO
514 END DO
515 END DO
516 CALL timestop(handle_sf)
517
518 !loop over angular sym fnct grps
519 CALL timeset('nnp_acsf_angular', handle_sf)
520 off = nnp%n_rad(ind)
521 DO s = 1, nnp%ang(ind)%n_symfgrp
522 cutoff_s = nnp%ang(ind)%symfgrp(s)%cutoff
523 cutoff_sqr = cutoff_s*cutoff_s
524 n_symf_s = nnp%ang(ind)%symfgrp(s)%n_symf
525 n_ang1_s = neighbor%n_ang1(s)
526
527 ! Precompute cutoff values for ang1 neighbors (no derivatives needed)
528 CALL nnp_fill_fc_cache(neighbor%ang1(s)%dist, n_ang1_s, &
529 nnp%cut_type, cutoff_s, fc_cache1)
530
531 IF (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2)) THEN
532 DO j = 1, n_ang1_s
533 rvect1 = neighbor%ang1(s)%dist(1:3, j)
534 r1 = neighbor%ang1(s)%dist(4, j)
535 DO k = j + 1, n_ang1_s
536 rvect2 = neighbor%ang1(s)%dist(1:3, k)
537 r2 = neighbor%ang1(s)%dist(4, k)
538 rvect3(1) = rvect2(1) - rvect1(1)
539 rvect3(2) = rvect2(2) - rvect1(2)
540 rvect3(3) = rvect2(3) - rvect1(3)
541 r3_sqr = rvect3(1)*rvect3(1) + rvect3(2)*rvect3(2) + rvect3(3)*rvect3(3)
542 IF (r3_sqr < cutoff_sqr) THEN
543 r3 = sqrt(r3_sqr)
544 CALL nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, r1, r2, r3, &
545 fc_cache1(j), 0.0_dp, fc_cache1(k), 0.0_dp, &
546 angular_symtmp(1:n_symf_s))
547 DO sf = 1, n_symf_s
548 m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
549 nnp%ang(ind)%y(m - off) = nnp%ang(ind)%y(m - off) + angular_symtmp(sf)
550 END DO
551 END IF
552 END DO
553 END DO
554 ELSE
555 ! Precompute cutoff values for ang2 neighbors
556 n_ang2_s = neighbor%n_ang2(s)
557 CALL nnp_fill_fc_cache(neighbor%ang2(s)%dist, n_ang2_s, &
558 nnp%cut_type, cutoff_s, fc_cache2)
559
560 DO j = 1, n_ang1_s
561 rvect1 = neighbor%ang1(s)%dist(1:3, j)
562 r1 = neighbor%ang1(s)%dist(4, j)
563 DO k = 1, n_ang2_s
564 rvect2 = neighbor%ang2(s)%dist(1:3, k)
565 r2 = neighbor%ang2(s)%dist(4, k)
566 rvect3(1) = rvect2(1) - rvect1(1)
567 rvect3(2) = rvect2(2) - rvect1(2)
568 rvect3(3) = rvect2(3) - rvect1(3)
569 r3_sqr = rvect3(1)*rvect3(1) + rvect3(2)*rvect3(2) + rvect3(3)*rvect3(3)
570 IF (r3_sqr < cutoff_sqr) THEN
571 r3 = sqrt(r3_sqr)
572 CALL nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, r1, r2, r3, &
573 fc_cache1(j), 0.0_dp, fc_cache2(k), 0.0_dp, &
574 angular_symtmp(1:n_symf_s))
575 DO sf = 1, n_symf_s
576 m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
577 nnp%ang(ind)%y(m - off) = nnp%ang(ind)%y(m - off) + angular_symtmp(sf)
578 END DO
579 END IF
580 END DO
581 END DO
582 END IF
583 END DO
584 CALL timestop(handle_sf)
585 END IF
586
587 END associate
588
589 ! fc_cache1/2 and dfc_cache1/2 are persistent workspace; nothing to deallocate here.
590
591 END associate
592
593 !check extrapolation
594 CALL nnp_check_extrapolation(nnp, ind)
595
596 IF (PRESENT(stress)) THEN
597 CALL nnp_scale_acsf(nnp, ind, do_forces, stress)
598 ELSE
599 CALL nnp_scale_acsf(nnp, ind, do_forces)
600 END IF
601
602 CALL timestop(handle)
603
604 END SUBROUTINE nnp_calc_acsf
605
606! **************************************************************************************************
607!> \brief Fill the per-neighbour fc and dfc cutoff-function caches used by the
608!> force branch of the angular ACSF kernel.
609!> \param dist (4, :) neighbour array; column 4 holds the scalar distance
610!> \param n number of neighbors to process
611!> \param cut_type cutoff function selector (nnp_cut_cos / nnp_cut_tanh)
612!> \param cutoff_s per-group cutoff radius
613!> \param fc_cache output fc values, sized >= n
614!> \param dfc_cache output dfc values, sized >= n
615! **************************************************************************************************
616 PURE SUBROUTINE nnp_fill_fc_dfc_cache(dist, n, cut_type, cutoff_s, fc_cache, dfc_cache)
617 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: dist
618 INTEGER, INTENT(IN) :: n, cut_type
619 REAL(kind=dp), INTENT(IN) :: cutoff_s
620 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: fc_cache, dfc_cache
621
622 INTEGER :: j
623 REAL(kind=dp) :: arg_tmp, r_tmp, tanh_tmp
624
625 DO j = 1, n
626 r_tmp = dist(4, j)
627 SELECT CASE (cut_type)
628 CASE (nnp_cut_cos)
629 arg_tmp = pi*r_tmp/cutoff_s
630 fc_cache(j) = 0.5_dp*(cos(arg_tmp) + 1.0_dp)
631 dfc_cache(j) = -0.5_dp*sin(arg_tmp)*(pi/cutoff_s)
632 CASE (nnp_cut_tanh)
633 tanh_tmp = tanh(1.0_dp - r_tmp/cutoff_s)
634 fc_cache(j) = tanh_tmp**3
635 dfc_cache(j) = (-3.0_dp/cutoff_s)*(tanh_tmp**2 - tanh_tmp**4)
636 END SELECT
637 END DO
638
639 END SUBROUTINE nnp_fill_fc_dfc_cache
640
641! **************************************************************************************************
642!> \brief Fill the per-neighbour fc cutoff-function cache for the sym-only
643!> (no-forces) branch of the angular ACSF kernel. Derivatives are not needed.
644!> \param dist (4, :) neighbour array; column 4 holds the scalar distance
645!> \param n number of neighbors to process
646!> \param cut_type cutoff function selector (nnp_cut_cos / nnp_cut_tanh)
647!> \param cutoff_s per-group cutoff radius
648!> \param fc_cache output fc values, sized >= n
649! **************************************************************************************************
650 PURE SUBROUTINE nnp_fill_fc_cache(dist, n, cut_type, cutoff_s, fc_cache)
651 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: dist
652 INTEGER, INTENT(IN) :: n, cut_type
653 REAL(kind=dp), INTENT(IN) :: cutoff_s
654 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: fc_cache
655
656 INTEGER :: j
657 REAL(kind=dp) :: r_tmp, tanh_tmp
658
659 DO j = 1, n
660 r_tmp = dist(4, j)
661 SELECT CASE (cut_type)
662 CASE (nnp_cut_cos)
663 fc_cache(j) = 0.5_dp*(cos(pi*r_tmp/cutoff_s) + 1.0_dp)
664 CASE (nnp_cut_tanh)
665 tanh_tmp = tanh(1.0_dp - r_tmp/cutoff_s)
666 fc_cache(j) = tanh_tmp**3
667 END SELECT
668 END DO
669
670 END SUBROUTINE nnp_fill_fc_cache
671
672! **************************************************************************************************
673!> \brief OpenMP parallelization of the angular SF group loop over s. Groups
674!> partition the angular SF indices disjointly, so per-group writes into
675!> self_dGdr / stress / nnp%ang%y and the workspace dGdr_ang accumulators
676!> are race-free. Per-thread scratch is PRIVATE automatic arrays; buffer
677!> growth and zero-init run in a serial pre-pass so the parallel region
678!> never touches ALLOCATABLE state.
679!> \param nnp ...
680!> \param ind ...
681!> \param self_dGdr ...
682!> \param off ...
683!> \param stress ...
684! **************************************************************************************************
685 SUBROUTINE nnp_acsf_angular_loop_omp(nnp, ind, self_dGdr, off, stress)
686 TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
687 INTEGER, INTENT(IN) :: ind
688 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: self_dgdr
689 INTEGER, INTENT(IN) :: off
690 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
691 OPTIONAL :: stress
692
693 INTEGER :: cache_cap_loc, j, k, l, m, &
694 max_ang_symf_loc, n_ang1_s, n_ang2_s, &
695 n_symf_s, s, sf
696 LOGICAL :: homo_grp
697 REAL(kind=dp) :: cutoff_s, cutoff_sqr, r1, r2, r3, r3_sqr
698 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dfc_c1_loc, dfc_c2_loc, fc_c1_loc, &
699 fc_c2_loc, sym_loc
700 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: force_loc
701 REAL(kind=dp), DIMENSION(3) :: rvect1, rvect2, rvect3
702
703! Per-thread automatic scratch (see PRIVATE clause below). Sizes
704! are pulled from the pre-sized workspace so they match the peak
705! usage that the serial kernel would allocate into the workspace
706! %fc_cache*/angular_* slabs.
707
708 associate(workspace => nnp%neighbor_interface_state%workspace(ind), &
709 neighbor => nnp%neighbor_interface_state%workspace(ind)%neighbor)
710
711 cache_cap_loc = max(1, workspace%cache_cap)
712 max_ang_symf_loc = max(1, workspace%max_ang_symf)
713
714 ! Serial pre-pass: grow per-group accumulator slabs and zero them
715 ! out. Doing this outside the parallel region guarantees no thread
716 ! ever touches ALLOCATABLE components of nnp_dGdr_grp_type.
717 DO s = 1, nnp%ang(ind)%n_symfgrp
718 n_symf_s = nnp%ang(ind)%symfgrp(s)%n_symf
719 n_ang1_s = neighbor%n_ang1(s)
720 homo_grp = (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2))
721 CALL nnp_grp_grow_dgdr(workspace%dGdr_ang_jj(s), n_ang1_s)
722 IF (homo_grp) THEN
723 CALL nnp_grp_grow_dgdr(workspace%dGdr_ang_kk(s), n_ang1_s)
724 n_ang2_s = 0
725 ELSE
726 n_ang2_s = neighbor%n_ang2(s)
727 CALL nnp_grp_grow_dgdr(workspace%dGdr_ang_kk(s), n_ang2_s)
728 END IF
729 IF (n_ang1_s > 0) workspace%dGdr_ang_jj(s)%data(:, 1:n_symf_s, 1:n_ang1_s) = 0.0_dp
730 IF (homo_grp) THEN
731 IF (n_ang1_s > 0) workspace%dGdr_ang_kk(s)%data(:, 1:n_symf_s, 1:n_ang1_s) = 0.0_dp
732 ELSE
733 IF (n_ang2_s > 0) workspace%dGdr_ang_kk(s)%data(:, 1:n_symf_s, 1:n_ang2_s) = 0.0_dp
734 END IF
735 END DO
736
737 ! Each thread allocates its own PRIVATE scratch inside the parallel
738 ! region: an ALLOCATABLE listed as PRIVATE enters unallocated per
739 ! thread, so the ALLOCATE below gives each thread an independent slab.
740 ! Angular groups partition the SF index disjointly, so the shared
741 ! accumulators (self_dGdr, stress, nnp%ang%y, workspace dGdr) are written
742 ! race-free. workspace/neighbor are ASSOCIATE names that inherit the
743 ! data-sharing of their nnp selector; the OPTIONAL stress is SHARED and
744 ! gated by PRESENT().
745!$OMP PARALLEL DEFAULT(NONE) &
746!$OMP SHARED(nnp, ind, self_dGdr, off, stress, &
747!$OMP cache_cap_loc, max_ang_symf_loc) &
748!$OMP PRIVATE(s, j, k, sf, m, l, n_symf_s, n_ang1_s, n_ang2_s, &
749!$OMP r1, r2, r3, r3_sqr, cutoff_s, cutoff_sqr, homo_grp, &
750!$OMP rvect1, rvect2, rvect3, &
751!$OMP fc_c1_loc, dfc_c1_loc, fc_c2_loc, dfc_c2_loc, &
752!$OMP sym_loc, force_loc)
753 ALLOCATE (fc_c1_loc(cache_cap_loc))
754 ALLOCATE (dfc_c1_loc(cache_cap_loc))
755 ALLOCATE (fc_c2_loc(cache_cap_loc))
756 ALLOCATE (dfc_c2_loc(cache_cap_loc))
757 ALLOCATE (sym_loc(max_ang_symf_loc))
758 ALLOCATE (force_loc(3, 3, max_ang_symf_loc))
759
760!$OMP DO SCHEDULE(dynamic)
761 DO s = 1, nnp%ang(ind)%n_symfgrp
762 cutoff_s = nnp%ang(ind)%symfgrp(s)%cutoff
763 cutoff_sqr = cutoff_s*cutoff_s
764 n_symf_s = nnp%ang(ind)%symfgrp(s)%n_symf
765 n_ang1_s = neighbor%n_ang1(s)
766 homo_grp = (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2))
767
768 CALL nnp_fill_fc_dfc_cache(neighbor%ang1(s)%dist, n_ang1_s, &
769 nnp%cut_type, cutoff_s, fc_c1_loc, dfc_c1_loc)
770
771 IF (homo_grp) THEN
772 DO j = 1, n_ang1_s
773 rvect1 = neighbor%ang1(s)%dist(1:3, j)
774 r1 = neighbor%ang1(s)%dist(4, j)
775 DO k = j + 1, n_ang1_s
776 rvect2 = neighbor%ang1(s)%dist(1:3, k)
777 r2 = neighbor%ang1(s)%dist(4, k)
778 rvect3(1) = rvect2(1) - rvect1(1)
779 rvect3(2) = rvect2(2) - rvect1(2)
780 rvect3(3) = rvect2(3) - rvect1(3)
781 r3_sqr = rvect3(1)*rvect3(1) + rvect3(2)*rvect3(2) + rvect3(3)*rvect3(3)
782 IF (r3_sqr < cutoff_sqr) THEN
783 r3 = sqrt(r3_sqr)
784 CALL nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, &
785 r1, r2, r3, &
786 fc_c1_loc(j), dfc_c1_loc(j), &
787 fc_c1_loc(k), dfc_c1_loc(k), &
788 sym_loc(1:n_symf_s), &
789 force_loc(:, :, 1:n_symf_s))
790 DO sf = 1, n_symf_s
791 m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
792 self_dgdr(1, m) = self_dgdr(1, m) + force_loc(1, 1, sf)
793 self_dgdr(2, m) = self_dgdr(2, m) + force_loc(2, 1, sf)
794 self_dgdr(3, m) = self_dgdr(3, m) + force_loc(3, 1, sf)
795 workspace%dGdr_ang_jj(s)%data(1, sf, j) = workspace%dGdr_ang_jj(s)%data(1, sf, j) + force_loc(1, 2, sf)
796 workspace%dGdr_ang_jj(s)%data(2, sf, j) = workspace%dGdr_ang_jj(s)%data(2, sf, j) + force_loc(2, 2, sf)
797 workspace%dGdr_ang_jj(s)%data(3, sf, j) = workspace%dGdr_ang_jj(s)%data(3, sf, j) + force_loc(3, 2, sf)
798 workspace%dGdr_ang_kk(s)%data(1, sf, k) = workspace%dGdr_ang_kk(s)%data(1, sf, k) + force_loc(1, 3, sf)
799 workspace%dGdr_ang_kk(s)%data(2, sf, k) = workspace%dGdr_ang_kk(s)%data(2, sf, k) + force_loc(2, 3, sf)
800 workspace%dGdr_ang_kk(s)%data(3, sf, k) = workspace%dGdr_ang_kk(s)%data(3, sf, k) + force_loc(3, 3, sf)
801 IF (PRESENT(stress)) THEN
802 DO l = 1, 3
803 stress(:, l, m) = stress(:, l, m) - rvect1(:)*force_loc(l, 2, sf)
804 stress(:, l, m) = stress(:, l, m) - rvect2(:)*force_loc(l, 3, sf)
805 END DO
806 END IF
807 nnp%ang(ind)%y(m - off) = nnp%ang(ind)%y(m - off) + sym_loc(sf)
808 END DO
809 END IF
810 END DO
811 END DO
812 ELSE
813 n_ang2_s = neighbor%n_ang2(s)
814 CALL nnp_fill_fc_dfc_cache(neighbor%ang2(s)%dist, n_ang2_s, &
815 nnp%cut_type, cutoff_s, fc_c2_loc, dfc_c2_loc)
816
817 DO j = 1, n_ang1_s
818 rvect1 = neighbor%ang1(s)%dist(1:3, j)
819 r1 = neighbor%ang1(s)%dist(4, j)
820 DO k = 1, n_ang2_s
821 rvect2 = neighbor%ang2(s)%dist(1:3, k)
822 r2 = neighbor%ang2(s)%dist(4, k)
823 rvect3(1) = rvect2(1) - rvect1(1)
824 rvect3(2) = rvect2(2) - rvect1(2)
825 rvect3(3) = rvect2(3) - rvect1(3)
826 r3_sqr = rvect3(1)*rvect3(1) + rvect3(2)*rvect3(2) + rvect3(3)*rvect3(3)
827 IF (r3_sqr < cutoff_sqr) THEN
828 r3 = sqrt(r3_sqr)
829 CALL nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, &
830 r1, r2, r3, &
831 fc_c1_loc(j), dfc_c1_loc(j), &
832 fc_c2_loc(k), dfc_c2_loc(k), &
833 sym_loc(1:n_symf_s), &
834 force_loc(:, :, 1:n_symf_s))
835 DO sf = 1, n_symf_s
836 m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
837 self_dgdr(1, m) = self_dgdr(1, m) + force_loc(1, 1, sf)
838 self_dgdr(2, m) = self_dgdr(2, m) + force_loc(2, 1, sf)
839 self_dgdr(3, m) = self_dgdr(3, m) + force_loc(3, 1, sf)
840 workspace%dGdr_ang_jj(s)%data(1, sf, j) = workspace%dGdr_ang_jj(s)%data(1, sf, j) + force_loc(1, 2, sf)
841 workspace%dGdr_ang_jj(s)%data(2, sf, j) = workspace%dGdr_ang_jj(s)%data(2, sf, j) + force_loc(2, 2, sf)
842 workspace%dGdr_ang_jj(s)%data(3, sf, j) = workspace%dGdr_ang_jj(s)%data(3, sf, j) + force_loc(3, 2, sf)
843 workspace%dGdr_ang_kk(s)%data(1, sf, k) = workspace%dGdr_ang_kk(s)%data(1, sf, k) + force_loc(1, 3, sf)
844 workspace%dGdr_ang_kk(s)%data(2, sf, k) = workspace%dGdr_ang_kk(s)%data(2, sf, k) + force_loc(2, 3, sf)
845 workspace%dGdr_ang_kk(s)%data(3, sf, k) = workspace%dGdr_ang_kk(s)%data(3, sf, k) + force_loc(3, 3, sf)
846 IF (PRESENT(stress)) THEN
847 DO l = 1, 3
848 stress(:, l, m) = stress(:, l, m) - rvect1(:)*force_loc(l, 2, sf)
849 stress(:, l, m) = stress(:, l, m) - rvect2(:)*force_loc(l, 3, sf)
850 END DO
851 END IF
852 nnp%ang(ind)%y(m - off) = nnp%ang(ind)%y(m - off) + sym_loc(sf)
853 END DO
854 END IF
855 END DO
856 END DO
857 END IF
858 END DO
859!$OMP END DO
860
861 DEALLOCATE (fc_c1_loc, dfc_c1_loc, fc_c2_loc, dfc_c2_loc, sym_loc, force_loc)
862!$OMP END PARALLEL
863
864 END associate
865
866 END SUBROUTINE nnp_acsf_angular_loop_omp
867
868! **************************************************************************************************
869!> \brief Prepare or update the linked-cell / Verlet cache for the current
870!> geometry. Lazily allocates nnp%cell_list_cache and
871!> nnp%neighbor_interface_state, then delegates to
872!> nnp_prepare_cell_list_cache and nnp_neighbor_interface_prepare. Call
873!> once per force eval before the per-atom loop; cheap on re-entry when
874!> nothing has changed.
875!> \param nnp NNP environment whose persistent caches are to be (re)built.
876!> \author Dhruv Sharma (ds2173@cam.ac.uk)
877! **************************************************************************************************
878 SUBROUTINE nnp_prepare_neighbor_cache(nnp)
879
880 TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
881
882 IF (.NOT. ALLOCATED(nnp%cell_list_cache)) ALLOCATE (nnp%cell_list_cache)
883 IF (.NOT. ALLOCATED(nnp%neighbor_interface_state)) ALLOCATE (nnp%neighbor_interface_state)
884 CALL nnp_prepare_cell_list_cache(nnp)
885 CALL nnp_neighbor_interface_prepare(nnp)
886
887 END SUBROUTINE nnp_prepare_neighbor_cache
888
889! **************************************************************************************************
890!> \brief Check if the nnp is extrapolating
891!> \param nnp ...
892!> \param ind ...
893!> \date 2020-10-10
894!> \author Christoph Schran (christoph.schran@rub.de)
895! **************************************************************************************************
896 SUBROUTINE nnp_check_extrapolation(nnp, ind)
897
898 TYPE(nnp_type), INTENT(INOUT) :: nnp
899 INTEGER, INTENT(IN) :: ind
900
901 REAL(kind=dp), PARAMETER :: threshold = 0.0001_dp
902
903 INTEGER :: j
904 LOGICAL :: extrapolate
905
906 extrapolate = nnp%output_expol
907
908 DO j = 1, nnp%n_rad(ind)
909 IF (nnp%rad(ind)%y(j) - nnp%rad(ind)%loc_max(j) > threshold) THEN
910 extrapolate = .true.
911 ELSE IF (-nnp%rad(ind)%y(j) + nnp%rad(ind)%loc_min(j) > threshold) THEN
912 extrapolate = .true.
913 END IF
914 END DO
915 DO j = 1, nnp%n_ang(ind)
916 IF (nnp%ang(ind)%y(j) - nnp%ang(ind)%loc_max(j) > threshold) THEN
917 extrapolate = .true.
918 ELSE IF (-nnp%ang(ind)%y(j) + nnp%ang(ind)%loc_min(j) > threshold) THEN
919 extrapolate = .true.
920 END IF
921 END DO
922
923 nnp%output_expol = extrapolate
924
925 END SUBROUTINE nnp_check_extrapolation
926
927! **************************************************************************************************
928!> \brief Scale and center symmetry functions (and gradients)
929!> \param nnp ...
930!> \param ind ...
931!> \param do_forces ...
932!> \param stress ...
933!> \date 2020-10-10
934!> \author Christoph Schran (christoph.schran@rub.de)
935! **************************************************************************************************
936 SUBROUTINE nnp_scale_acsf(nnp, ind, do_forces, stress)
937
938 TYPE(nnp_type), INTENT(INOUT) :: nnp
939 INTEGER, INTENT(IN) :: ind
940 LOGICAL, INTENT(IN) :: do_forces
941 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
942 OPTIONAL :: stress
943
944 INTEGER :: j, k, m, n_ang1_s, n_ang2_s, n_symf_s, &
945 off, s, sf
946 LOGICAL :: homo_grp
947 REAL(kind=dp) :: scale
948
949! INOUT (not OUT): stress is per-input-node and accumulates across central atoms.
950
951 IF (nnp%center_acsf) THEN
952 DO j = 1, nnp%n_rad(ind)
953 nnp%arc(ind)%layer(1)%node(j) = nnp%rad(ind)%y(j) - nnp%rad(ind)%loc_av(j)
954 END DO
955 off = nnp%n_rad(ind)
956 DO j = 1, nnp%n_ang(ind)
957 nnp%arc(ind)%layer(1)%node(j + off) = nnp%ang(ind)%y(j) - nnp%ang(ind)%loc_av(j)
958 END DO
959
960 IF (nnp%scale_acsf) THEN
961 DO j = 1, nnp%n_rad(ind)
962 nnp%arc(ind)%layer(1)%node(j) = nnp%arc(ind)%layer(1)%node(j)/ &
963 (nnp%rad(ind)%loc_max(j) - nnp%rad(ind)%loc_min(j))*(nnp%scmax - nnp%scmin) + nnp%scmin
964 END DO
965 off = nnp%n_rad(ind)
966 DO j = 1, nnp%n_ang(ind)
967 nnp%arc(ind)%layer(1)%node(j + off) = nnp%arc(ind)%layer(1)%node(j + off)/ &
968 (nnp%ang(ind)%loc_max(j) - nnp%ang(ind)%loc_min(j))*(nnp%scmax - nnp%scmin) + nnp%scmin
969 END DO
970 END IF
971 ELSE IF (nnp%scale_acsf) THEN
972 DO j = 1, nnp%n_rad(ind)
973 nnp%arc(ind)%layer(1)%node(j) = (nnp%rad(ind)%y(j) - nnp%rad(ind)%loc_min(j))/ &
974 (nnp%rad(ind)%loc_max(j) - nnp%rad(ind)%loc_min(j))* &
975 (nnp%scmax - nnp%scmin) + nnp%scmin
976 END DO
977 off = nnp%n_rad(ind)
978 DO j = 1, nnp%n_ang(ind)
979 nnp%arc(ind)%layer(1)%node(j + off) = (nnp%ang(ind)%y(j) - nnp%ang(ind)%loc_min(j))/ &
980 (nnp%ang(ind)%loc_max(j) - nnp%ang(ind)%loc_min(j))* &
981 (nnp%scmax - nnp%scmin) + nnp%scmin
982 END DO
983 ELSE IF (nnp%scale_sigma_acsf) THEN
984 DO j = 1, nnp%n_rad(ind)
985 nnp%arc(ind)%layer(1)%node(j) = (nnp%rad(ind)%y(j) - nnp%rad(ind)%loc_av(j))/ &
986 nnp%rad(ind)%sigma(j)*(nnp%scmax - nnp%scmin) + nnp%scmin
987 END DO
988 off = nnp%n_rad(ind)
989 DO j = 1, nnp%n_ang(ind)
990 nnp%arc(ind)%layer(1)%node(j + off) = (nnp%ang(ind)%y(j) - nnp%ang(ind)%loc_av(j))/ &
991 nnp%ang(ind)%sigma(j)*(nnp%scmax - nnp%scmin) + nnp%scmin
992 END DO
993 ELSE
994 DO j = 1, nnp%n_rad(ind)
995 nnp%arc(ind)%layer(1)%node(j) = nnp%rad(ind)%y(j)
996 END DO
997 off = nnp%n_rad(ind)
998 DO j = 1, nnp%n_ang(ind)
999 nnp%arc(ind)%layer(1)%node(j + off) = nnp%ang(ind)%y(j)
1000 END DO
1001 END IF
1002
1003 IF (do_forces .AND. (nnp%scale_acsf .OR. nnp%scale_sigma_acsf)) THEN
1004 ! Scale the per-neighbor dGdr slabs in workspace state. Only the
1005 ! actually-populated entries in each slab need touching: for each
1006 ! group s the valid ranges are sf=1..n_symf_s and j=1..n_rad(s)
1007 ! (radial) or j=1..n_ang1_s, k=1..n_ang2_s (angular). The self
1008 ! contribution lives in self_dGdr(:, m).
1009 associate(workspace => nnp%neighbor_interface_state%workspace(ind), &
1010 neighbor => nnp%neighbor_interface_state%workspace(ind)%neighbor, &
1011 self_dgdr => nnp%neighbor_interface_state%workspace(ind)%self_dGdr)
1012
1013 ! Radial groups
1014 DO s = 1, nnp%rad(ind)%n_symfgrp
1015 n_symf_s = nnp%rad(ind)%symfgrp(s)%n_symf
1016 associate(rad_buf => workspace%dGdr_rad(s)%data)
1017 DO sf = 1, n_symf_s
1018 m = nnp%rad(ind)%symfgrp(s)%symf(sf)
1019 IF (nnp%scale_acsf) THEN
1020 scale = (nnp%scmax - nnp%scmin)/ &
1021 (nnp%rad(ind)%loc_max(m) - nnp%rad(ind)%loc_min(m))
1022 ELSE
1023 scale = (nnp%scmax - nnp%scmin)/nnp%rad(ind)%sigma(m)
1024 END IF
1025 self_dgdr(1, m) = self_dgdr(1, m)*scale
1026 self_dgdr(2, m) = self_dgdr(2, m)*scale
1027 self_dgdr(3, m) = self_dgdr(3, m)*scale
1028 DO j = 1, neighbor%n_rad(s)
1029 rad_buf(1, sf, j) = rad_buf(1, sf, j)*scale
1030 rad_buf(2, sf, j) = rad_buf(2, sf, j)*scale
1031 rad_buf(3, sf, j) = rad_buf(3, sf, j)*scale
1032 END DO
1033 END DO
1034 END associate
1035 END DO
1036
1037 ! Angular groups
1038 off = nnp%n_rad(ind)
1039 DO s = 1, nnp%ang(ind)%n_symfgrp
1040 n_symf_s = nnp%ang(ind)%symfgrp(s)%n_symf
1041 n_ang1_s = neighbor%n_ang1(s)
1042 homo_grp = (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2))
1043 IF (homo_grp) THEN
1044 ! kk slab is also indexed in ind_ang1 in the homo case
1045 n_ang2_s = n_ang1_s
1046 ELSE
1047 n_ang2_s = neighbor%n_ang2(s)
1048 END IF
1049 associate(jj_buf => workspace%dGdr_ang_jj(s)%data, &
1050 kk_buf => workspace%dGdr_ang_kk(s)%data)
1051 DO sf = 1, n_symf_s
1052 m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
1053 IF (nnp%scale_acsf) THEN
1054 scale = (nnp%scmax - nnp%scmin)/ &
1055 (nnp%ang(ind)%loc_max(m - off) - nnp%ang(ind)%loc_min(m - off))
1056 ELSE
1057 scale = (nnp%scmax - nnp%scmin)/nnp%ang(ind)%sigma(m - off)
1058 END IF
1059 self_dgdr(1, m) = self_dgdr(1, m)*scale
1060 self_dgdr(2, m) = self_dgdr(2, m)*scale
1061 self_dgdr(3, m) = self_dgdr(3, m)*scale
1062 DO j = 1, n_ang1_s
1063 jj_buf(1, sf, j) = jj_buf(1, sf, j)*scale
1064 jj_buf(2, sf, j) = jj_buf(2, sf, j)*scale
1065 jj_buf(3, sf, j) = jj_buf(3, sf, j)*scale
1066 END DO
1067 DO k = 1, n_ang2_s
1068 kk_buf(1, sf, k) = kk_buf(1, sf, k)*scale
1069 kk_buf(2, sf, k) = kk_buf(2, sf, k)*scale
1070 kk_buf(3, sf, k) = kk_buf(3, sf, k)*scale
1071 END DO
1072 END DO
1073 END associate
1074 END DO
1075
1076 END associate
1077 END IF
1078
1079 IF (PRESENT(stress)) THEN
1080 IF (nnp%scale_acsf) THEN
1081 DO j = 1, nnp%n_rad(ind)
1082 stress(:, :, j) = stress(:, :, j)/(nnp%rad(ind)%loc_max(j) - nnp%rad(ind)%loc_min(j))* &
1083 (nnp%scmax - nnp%scmin)
1084 END DO
1085 off = nnp%n_rad(ind)
1086 DO j = 1, nnp%n_ang(ind)
1087 stress(:, :, j + off) = stress(:, :, j + off)/ &
1088 (nnp%ang(ind)%loc_max(j) - nnp%ang(ind)%loc_min(j))* &
1089 (nnp%scmax - nnp%scmin)
1090 END DO
1091 ELSE IF (nnp%scale_sigma_acsf) THEN
1092 DO j = 1, nnp%n_rad(ind)
1093 stress(:, :, j) = stress(:, :, j)/nnp%rad(ind)%sigma(j)*(nnp%scmax - nnp%scmin)
1094 END DO
1095 off = nnp%n_rad(ind)
1096 DO j = 1, nnp%n_ang(ind)
1097 stress(:, :, j + off) = stress(:, :, j + off)/nnp%ang(ind)%sigma(j)*(nnp%scmax - nnp%scmin)
1098 END DO
1099 END IF
1100 END IF
1101
1102 END SUBROUTINE nnp_scale_acsf
1103
1104! **************************************************************************************************
1105!> \brief Calculate radial symmetry function and gradient (optional)
1106!> \param nnp ...
1107!> \param ind ...
1108!> \param s ...
1109!> \param rvect ...
1110!> \param r ...
1111!> \param sym ...
1112!> \param force ...
1113!> \date 2020-10-10
1114!> \author Christoph Schran (christoph.schran@rub.de)
1115! **************************************************************************************************
1116 SUBROUTINE nnp_calc_rad(nnp, ind, s, rvect, r, sym, force)
1117
1118 TYPE(nnp_type), INTENT(IN), TARGET :: nnp
1119 INTEGER, INTENT(IN) :: ind, s
1120 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rvect
1121 REAL(kind=dp), INTENT(IN) :: r
1122 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: sym
1123 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT), &
1124 OPTIONAL :: force
1125
1126 INTEGER :: i, n_symf, sf
1127 REAL(kind=dp) :: dh00, dh01, dh10, dh11, drdx_x, drdx_y, &
1128 drdx_z, dsymdr_full, dyi, dyi1, h00, &
1129 h01, h10, h10_dx, h11, h11_dx, r_inv, &
1130 t, t2, t3, yi, yi1
1131 TYPE(nnp_symfgrp_type), POINTER :: grp
1132
1133 grp => nnp%rad(ind)%symfgrp(s)
1134 n_symf = grp%n_symf
1135
1136 ! Group-shared Hermite cubic spline. All SFs in this radial group share
1137 ! grp%cutoff, hence the same uniform grid, so the interpolation parameters
1138 ! (i, t, h00..h11, dh00..dh11) are computed once outside the SF loop; the
1139 ! inner loop then streams 4 contiguous reads per SF from spline_y/spline_dy.
1140 !
1141 ! Out-of-range clamping: for r past spline_x_max (= grp%cutoff) every SF
1142 ! sees the boundary value of its tabulated y(r), which the build routine
1143 ! pinned to 0 since fcut(cutoff) = 0, so the clamp returns sym = 0.
1144 IF (r >= grp%spline_x_max) THEN
1145 DO sf = 1, n_symf
1146 sym(sf) = 0.0_dp
1147 END DO
1148 IF (PRESENT(force)) THEN
1149 DO sf = 1, n_symf
1150 force(1, sf) = 0.0_dp
1151 force(2, sf) = 0.0_dp
1152 force(3, sf) = 0.0_dp
1153 END DO
1154 END IF
1155 RETURN
1156 END IF
1157
1158 i = int(r*grp%spline_dx_inv) + 1
1159 IF (i < 1) i = 1
1160 IF (i > grp%spline_n - 1) i = grp%spline_n - 1
1161
1162 t = (r - real(i - 1, kind=dp)*grp%spline_dx)*grp%spline_dx_inv
1163 t2 = t*t
1164 t3 = t2*t
1165
1166 h00 = 2.0_dp*t3 - 3.0_dp*t2 + 1.0_dp
1167 h10 = t3 - 2.0_dp*t2 + t
1168 h01 = -2.0_dp*t3 + 3.0_dp*t2
1169 h11 = t3 - t2
1170 h10_dx = h10*grp%spline_dx
1171 h11_dx = h11*grp%spline_dx
1172
1173 IF (PRESENT(force)) THEN
1174 dh00 = 6.0_dp*(t2 - t)
1175 dh10 = 3.0_dp*t2 - 4.0_dp*t + 1.0_dp
1176 dh01 = -dh00
1177 dh11 = 3.0_dp*t2 - 2.0_dp*t
1178
1179 r_inv = 1.0_dp/r
1180 drdx_x = rvect(1)*r_inv
1181 drdx_y = rvect(2)*r_inv
1182 drdx_z = rvect(3)*r_inv
1183
1184 ! Vectorizable inner loop: contiguous reads on sf, no branches,
1185 ! no function calls, no transcendentals.
1186 associate(spy_i => grp%spline_y(:, i), spy_i1 => grp%spline_y(:, i + 1), &
1187 spdy_i => grp%spline_dy(:, i), spdy_i1 => grp%spline_dy(:, i + 1))
1188 !$OMP SIMD PRIVATE(yi, yi1, dyi, dyi1, dsymdr_full)
1189 DO sf = 1, n_symf
1190 yi = spy_i(sf)
1191 yi1 = spy_i1(sf)
1192 dyi = spdy_i(sf)
1193 dyi1 = spdy_i1(sf)
1194 sym(sf) = h00*yi + h10_dx*dyi + h01*yi1 + h11_dx*dyi1
1195 dsymdr_full = (dh00*yi + dh01*yi1)*grp%spline_dx_inv + dh10*dyi + dh11*dyi1
1196 force(1, sf) = dsymdr_full*drdx_x
1197 force(2, sf) = dsymdr_full*drdx_y
1198 force(3, sf) = dsymdr_full*drdx_z
1199 END DO
1200 END associate
1201 ELSE
1202 associate(spy_i => grp%spline_y(:, i), spy_i1 => grp%spline_y(:, i + 1), &
1203 spdy_i => grp%spline_dy(:, i), spdy_i1 => grp%spline_dy(:, i + 1))
1204 !$OMP SIMD
1205 DO sf = 1, n_symf
1206 sym(sf) = h00*spy_i(sf) + h10_dx*spdy_i(sf) + &
1207 h01*spy_i1(sf) + h11_dx*spdy_i1(sf)
1208 END DO
1209 END associate
1210 END IF
1211
1212 END SUBROUTINE nnp_calc_rad
1213
1214! **************************************************************************************************
1215!> \brief Build sf-first (n_symf, n_grid) Hermite cubic spline tables for radial SFs.
1216!> Uses a custom sf-first layout (not splines_methods) so nnp_calc_rad can
1217!> stream contiguous SF values under !$OMP SIMD without indirect addressing.
1218!> \param nnp ...
1219! **************************************************************************************************
1220 SUBROUTINE nnp_build_radial_splines(nnp)
1221
1222 TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
1223
1224 CHARACTER(len=*), PARAMETER :: routinen = 'nnp_build_radial_splines'
1225
1226 INTEGER :: handle, ind, k, n_symf, p, s, sf
1227 REAL(kind=dp) :: arg, cutoff, dfcutdr, dr, eta, exp_term, &
1228 fcut, r, rs, tanh_tmp
1229
1230 CALL timeset(routinen, handle)
1231
1232 DO ind = 1, nnp%n_ele
1233 DO s = 1, nnp%rad(ind)%n_symfgrp
1234 associate(grp => nnp%rad(ind)%symfgrp(s))
1235 n_symf = grp%n_symf
1236 cutoff = grp%cutoff
1237 dr = cutoff/real(nnp%rad_spline_n - 1, kind=dp)
1238
1239 IF (ALLOCATED(grp%spline_y)) DEALLOCATE (grp%spline_y)
1240 IF (ALLOCATED(grp%spline_dy)) DEALLOCATE (grp%spline_dy)
1241 ALLOCATE (grp%spline_y(max(1, n_symf), nnp%rad_spline_n))
1242 ALLOCATE (grp%spline_dy(max(1, n_symf), nnp%rad_spline_n))
1243 grp%spline_n = nnp%rad_spline_n
1244 grp%spline_dx = dr
1245 grp%spline_dx_inv = 1.0_dp/dr
1246 grp%spline_x_max = cutoff
1247
1248 ! Fill the table SF-by-SF. The grid loop is the inner one
1249 ! during build only; runtime nnp_calc_rad reads sf-first.
1250 DO sf = 1, n_symf
1251 k = grp%symf(sf)
1252 eta = nnp%rad(ind)%eta(k)
1253 rs = nnp%rad(ind)%rs(k)
1254
1255 DO p = 1, nnp%rad_spline_n
1256 r = real(p - 1, kind=dp)*dr
1257
1258 SELECT CASE (nnp%cut_type)
1259 CASE (nnp_cut_cos)
1260 arg = pi*r/cutoff
1261 fcut = 0.5_dp*(cos(arg) + 1.0_dp)
1262 dfcutdr = -0.5_dp*sin(arg)*(pi/cutoff)
1263 CASE (nnp_cut_tanh)
1264 tanh_tmp = tanh(1.0_dp - r/cutoff)
1265 fcut = tanh_tmp**3
1266 dfcutdr = (-3.0_dp/cutoff)*(tanh_tmp**2 - tanh_tmp**4)
1267 CASE DEFAULT
1268 cpabort("NNP| Cutoff function unknown")
1269 END SELECT
1270
1271 exp_term = exp(-eta*(r - rs)**2)
1272
1273 grp%spline_y(sf, p) = exp_term*fcut
1274 grp%spline_dy(sf, p) = exp_term*(-2.0_dp*eta*(r - rs))*fcut + &
1275 exp_term*dfcutdr
1276 END DO
1277 END DO
1278
1279 ! Pin the boundary so the runtime out-of-range branch can
1280 ! return zeros without re-checking each SF.
1281 DO sf = 1, n_symf
1282 grp%spline_y(sf, nnp%rad_spline_n) = 0.0_dp
1283 grp%spline_dy(sf, nnp%rad_spline_n) = 0.0_dp
1284 END DO
1285
1286 grp%spline_built = .true.
1287 END associate
1288 END DO
1289 END DO
1290
1291 CALL timestop(handle)
1292
1293 END SUBROUTINE nnp_build_radial_splines
1294
1295! **************************************************************************************************
1296!> \brief Calculate angular symmetry function and gradient (optional)
1297!>
1298!> Vectorized SF-batched form. The original eta-sorted exp-skip loop has
1299!> been replaced by a sequence of SF passes:
1300!> 1. Branchless angular base: clamps tmp at the cusp (tmp <= 0), then
1301!> computes tmpzeta via integer or real pow per SF. Stays scalar
1302!> because the int/real switch and integer pow cannot SIMD.
1303!> 2. Vectorized EXP: symtmp_arr(sf) = EXP(-eta(sf)*r2sum). Marked
1304!> !$OMP SIMD so the compiler maps it to libmvec/SVML vector EXP.
1305!> Vector libm hides EXP latency across the SF loop without needing
1306!> an eta-dedup pre-pass.
1307!> 3. SIMD sym scatter: sym(sf) = prefzeta*angular*symtmp*ftot.
1308!> 4. Force scatter (only when forces requested). The clamp in pass 1
1309!> forces tmpzeta_arr(sf) = 0 at the cusp, which makes both angular
1310!> AND pref_lam vanish naturally -- no per-iter branch.
1311!>
1312!> Caller passes precomputed cutoff values for the j/k legs; the j-k leg
1313!> is computed inline. Geometry-only derivative bases are hoisted outside
1314!> the SF loops as before.
1315!> \param nnp ...
1316!> \param ind ...
1317!> \param s ...
1318!> \param rvect1 ...
1319!> \param rvect2 ...
1320!> \param rvect3 ...
1321!> \param r1 ...
1322!> \param r2 ...
1323!> \param r3 ...
1324!> \param fcut_j ...
1325!> \param dfcut_j ...
1326!> \param fcut_k ...
1327!> \param dfcut_k ...
1328!> \param sym ...
1329!> \param force ...
1330!> \date 2020-10-10
1331!> \author Christoph Schran (christoph.schran@rub.de)
1332!> \note This kernel is also inlined into nnp_calc_acsf on the serial force
1333!> path; any change to the computation or force scatter here must be
1334!> mirrored at that inlining site (see the comment there).
1335! **************************************************************************************************
1336 SUBROUTINE nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, r1, r2, r3, &
1337 fcut_j, dfcut_j, fcut_k, dfcut_k, sym, force)
1338
1339 TYPE(nnp_type), INTENT(IN), TARGET :: nnp
1340 INTEGER, INTENT(IN) :: ind, s
1341 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rvect1, rvect2, rvect3
1342 REAL(kind=dp), INTENT(IN) :: r1, r2, r3, fcut_j, dfcut_j, fcut_k, &
1343 dfcut_k
1344 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: sym
1345 REAL(kind=dp), DIMENSION(:, :, :), INTENT(OUT), &
1346 OPTIONAL :: force
1347
1348 INTEGER :: ii, izeta, n_symf, sf
1349 LOGICAL :: do_forces
1350 REAL(kind=dp) :: angular, arg_tmp, costheta, dfcut3, dfcutdr1, dfcutdr2, dfcutdr3, dsymdr1, &
1351 dsymdr2, dsymdr3, eta, f, fcut3, fcut_rc, ftot, g, inv_g2, lam, pref_lam, prefzeta, &
1352 r2sum, rsqr1, rsqr2, rsqr3, symtmp, tanh_tmp, tmp, tmp1, tmp2, tmp3, tmpzeta, zeta
1353 REAL(kind=dp), DIMENSION(3) :: dcosbase1, dcosbase2, dcosbase3, dgdx1, &
1354 dgdx2, dgdx3, dr1dx, dr2dx, dr3dx
1355 TYPE(nnp_symfgrp_type), POINTER :: grp
1356 REAL(kind=dp), &
1357 DIMENSION(nnp%ang(ind)%symfgrp(s)%n_symf) :: angular_arr, symtmp_arr, tmpzeta_arr
1358
1359! Per-SF scratch (automatic arrays sized to the current group). For
1360! typical n2p2 ACSF n_symf is 8-32, so this is ~256-1024 bytes on
1361! stack per call. Holds the staged outputs of pass 1 + pass 2 so
1362! pass 3 streams them contiguously.
1363
1364 do_forces = PRESENT(force)
1365 grp => nnp%ang(ind)%symfgrp(s)
1366 n_symf = grp%n_symf
1367 fcut_rc = grp%cutoff
1368
1369 rsqr1 = r1*r1
1370 rsqr2 = r2*r2
1371 rsqr3 = r3*r3
1372 r2sum = rsqr1 + rsqr2 + rsqr3
1373
1374 f = rsqr3 - rsqr1 - rsqr2
1375 g = -2.0_dp*r1*r2
1376 costheta = f/g
1377
1378 ! Compute fcut3 for r3 (j-k distance -- cannot be precomputed by caller)
1379 SELECT CASE (nnp%cut_type)
1380 CASE (nnp_cut_cos)
1381 arg_tmp = pi*r3/fcut_rc
1382 fcut3 = 0.5_dp*(cos(arg_tmp) + 1.0_dp)
1383 IF (do_forces) dfcut3 = -0.5_dp*sin(arg_tmp)*(pi/fcut_rc)
1384 CASE (nnp_cut_tanh)
1385 tanh_tmp = tanh(1.0_dp - r3/fcut_rc)
1386 fcut3 = tanh_tmp**3
1387 IF (do_forces) dfcut3 = (-3.0_dp/fcut_rc)*(tanh_tmp**2 - tanh_tmp**4)
1388 CASE DEFAULT
1389 cpabort("NNP| Cutoff function unknown")
1390 END SELECT
1391
1392 ! Use precomputed fcut values for j and k neighbors
1393 ftot = fcut_j*fcut_k*fcut3
1394
1395 IF (do_forces) THEN
1396 ! Combined cutoff derivatives (product rule)
1397 dfcutdr1 = dfcut_j*fcut_k*fcut3
1398 dfcutdr2 = fcut_j*dfcut_k*fcut3
1399 dfcutdr3 = fcut_j*fcut_k*dfcut3
1400
1401 dr1dx(:) = rvect1(:)/r1
1402 dr2dx(:) = rvect2(:)/r2
1403 dr3dx(:) = rvect3(:)/r3
1404
1405 ! Hoist geometry-only parts of costheta derivatives outside SF loop.
1406 ! Full angular derivative: dangulardx = (zeta*tmpzeta*lam/g^2) * dcosbase
1407 ! where dcosbase factors out lam from dfdx and the f*dgdx term.
1408 ! inv_g2 = 1/g^2 is hoisted: one divide per triplet, one mul per SF.
1409 inv_g2 = 1.0_dp/(g*g)
1410 DO ii = 1, 3
1411 tmp1 = 2.0_dp*r2*dr1dx(ii)
1412 tmp2 = 2.0_dp*r1*dr2dx(ii)
1413 dgdx1(ii) = -(tmp1 + tmp2)
1414 dgdx2(ii) = tmp1
1415 dgdx3(ii) = tmp2
1416
1417 dcosbase1(ii) = -2.0_dp*(rvect1(ii) + rvect2(ii))*g - f*dgdx1(ii)
1418 dcosbase2(ii) = 2.0_dp*(rvect3(ii) + rvect1(ii))*g - f*dgdx2(ii)
1419 dcosbase3(ii) = 2.0_dp*(rvect2(ii) - rvect3(ii))*g - f*dgdx3(ii)
1420 END DO
1421 ELSE
1422 inv_g2 = 0.0_dp
1423 END IF
1424
1425 ! Pass 1: branchless cusp clamp -- zero tmpzeta propagates to sym and pref_lam.
1426 DO sf = 1, n_symf
1427 tmp = 1.0_dp + grp%pack_lam(sf)*costheta
1428 IF (tmp <= 0.0_dp) THEN
1429 tmpzeta_arr(sf) = 0.0_dp
1430 angular_arr(sf) = 0.0_dp
1431 ELSE
1432 IF (grp%pack_use_int_zeta(sf)) THEN
1433 izeta = grp%pack_izeta(sf)
1434 tmpzeta_arr(sf) = tmp**(izeta - 1)
1435 ELSE
1436 tmpzeta_arr(sf) = tmp**(grp%pack_zeta(sf) - 1.0_dp)
1437 END IF
1438 angular_arr(sf) = tmpzeta_arr(sf)*tmp
1439 END IF
1440 END DO
1441
1442 ! ---- Pass 2: vectorized EXP via libmvec/SVML ----
1443 ! Single SIMD loop computing symtmp_arr(sf) = EXP(-eta(sf)*r2sum).
1444 ! With -fopenmp-simd this maps to vector libm (e.g. _ZGVdN4v_exp on
1445 ! AVX2, _ZGVeN8v_exp on AVX-512), hiding EXP latency across 4-8 SFs
1446 ! per iter -- ~10 cycles/element vs ~30 for scalar EXP, which beats
1447 ! any eta-dedup pre-pass.
1448 !$OMP SIMD
1449 DO sf = 1, n_symf
1450 symtmp_arr(sf) = exp(-grp%pack_eta(sf)*r2sum)
1451 END DO
1452
1453 ! ---- Pass 3: sym scatter (vectorizable) ----
1454 !$OMP SIMD
1455 DO sf = 1, n_symf
1456 sym(sf) = grp%pack_prefzeta(sf)*angular_arr(sf)*symtmp_arr(sf)*ftot
1457 END DO
1458
1459 ! ---- Pass 4: force scatter (only when do_forces) ----
1460 IF (do_forces) THEN
1461 DO sf = 1, n_symf
1462 symtmp = symtmp_arr(sf)
1463 angular = angular_arr(sf)
1464 tmpzeta = tmpzeta_arr(sf)
1465 eta = grp%pack_eta(sf)
1466 lam = grp%pack_lam(sf)
1467 zeta = grp%pack_zeta(sf)
1468 prefzeta = grp%pack_prefzeta(sf)
1469
1470 ! pref_lam carries tmpzeta, so it vanishes at the cusp without
1471 ! an explicit branch (tmpzeta_arr was clamped to 0 in pass 1).
1472 pref_lam = zeta*tmpzeta*lam*inv_g2
1473
1474 tmp = -2.0_dp*symtmp*eta
1475 dsymdr1 = tmp*r1
1476 dsymdr2 = tmp*r2
1477 dsymdr3 = tmp*r3
1478
1479 tmp = prefzeta*symtmp*ftot
1480 tmp1 = prefzeta*angular*(ftot*dsymdr1 + dfcutdr1*symtmp)
1481 tmp2 = prefzeta*angular*(ftot*dsymdr2 + dfcutdr2*symtmp)
1482 tmp3 = prefzeta*angular*(ftot*dsymdr3 + dfcutdr3*symtmp)
1483 DO ii = 1, 3
1484 force(ii, 1, sf) = tmp*pref_lam*dcosbase1(ii) + tmp1*dr1dx(ii) + tmp2*dr2dx(ii)
1485 force(ii, 2, sf) = tmp*pref_lam*dcosbase2(ii) - tmp1*dr1dx(ii) + tmp3*dr3dx(ii)
1486 force(ii, 3, sf) = tmp*pref_lam*dcosbase3(ii) - tmp2*dr2dx(ii) - tmp3*dr3dx(ii)
1487 END DO
1488 END DO
1489 END IF
1490
1491 END SUBROUTINE nnp_calc_ang
1492
1493! **************************************************************************************************
1494!> \brief Sort an (ele, nuc_ele) pair of arrays in ascending order of atomic number.
1495!> Used to canonicalise element ordering inside the NNP environment so the
1496!> same model file produces the same per-element index layout regardless
1497!> of input ordering.
1498!> \param ele element-symbol array, sorted in place to match nuc_ele.
1499!> \param nuc_ele per-element atomic number, sorted in place.
1500!> \author Christoph Schran (christoph.schran@rub.de)
1501! **************************************************************************************************
1502 SUBROUTINE nnp_sort_ele(ele, nuc_ele)
1503 CHARACTER(len=2), DIMENSION(:), INTENT(INOUT) :: ele
1504 INTEGER, DIMENSION(:), INTENT(INOUT) :: nuc_ele
1505
1506 CHARACTER(len=2) :: tmp_ele
1507 INTEGER :: i, j, loc, minimum, tmp_nuc_ele
1508
1509 DO i = 1, SIZE(ele)
1510 CALL get_ptable_info(ele(i), number=nuc_ele(i))
1511 END DO
1512
1513 DO i = 1, SIZE(ele) - 1
1514 minimum = nuc_ele(i)
1515 loc = i
1516 DO j = i + 1, SIZE(ele)
1517 IF (nuc_ele(j) < minimum) THEN
1518 loc = j
1519 minimum = nuc_ele(j)
1520 END IF
1521 END DO
1522 tmp_nuc_ele = nuc_ele(i)
1523 nuc_ele(i) = nuc_ele(loc)
1524 nuc_ele(loc) = tmp_nuc_ele
1525
1526 tmp_ele = ele(i)
1527 ele(i) = ele(loc)
1528 ele(loc) = tmp_ele
1529 END DO
1530
1531 END SUBROUTINE nnp_sort_ele
1532
1533! **************************************************************************************************
1534!> \brief Sort radial and angular symmetry functions in canonical order.
1535!> Radial SFs are sorted by eta (ascending) then rcut; angular SFs by
1536!> eta, lambda, zeta. This is the order downstream code in nnp_init_acsf_groups
1537!> relies on for run-length-style group packing.
1538!> \param nnp NNP environment whose rad/ang SF arrays will be reordered in place.
1539!> \author Christoph Schran (christoph.schran@rub.de)
1540! **************************************************************************************************
1541 SUBROUTINE nnp_sort_acsf(nnp)
1542 TYPE(nnp_type), INTENT(INOUT) :: nnp
1543
1544 INTEGER :: i, j, k, loc
1545
1546 DO i = 1, nnp%n_ele
1547 DO j = 1, nnp%n_rad(i) - 1
1548 loc = j
1549 DO k = j + 1, nnp%n_rad(i)
1550 IF (nnp%rad(i)%funccut(loc) > nnp%rad(i)%funccut(k)) THEN
1551 loc = k
1552 END IF
1553 END DO
1554 CALL nnp_swaprad(nnp%rad(i), j, loc)
1555 END DO
1556
1557 DO j = 1, nnp%n_rad(i) - 1
1558 loc = j
1559 DO k = j + 1, nnp%n_rad(i)
1560 IF (nnp%rad(i)%funccut(loc) == nnp%rad(i)%funccut(k) .AND. &
1561 nnp%rad(i)%eta(loc) > nnp%rad(i)%eta(k)) THEN
1562 loc = k
1563 END IF
1564 END DO
1565 CALL nnp_swaprad(nnp%rad(i), j, loc)
1566 END DO
1567
1568 DO j = 1, nnp%n_rad(i) - 1
1569 loc = j
1570 DO k = j + 1, nnp%n_rad(i)
1571 IF (nnp%rad(i)%funccut(loc) == nnp%rad(i)%funccut(k) .AND. &
1572 nnp%rad(i)%eta(loc) == nnp%rad(i)%eta(k) .AND. &
1573 nnp%rad(i)%rs(loc) > nnp%rad(i)%rs(k)) THEN
1574 loc = k
1575 END IF
1576 END DO
1577 CALL nnp_swaprad(nnp%rad(i), j, loc)
1578 END DO
1579
1580 DO j = 1, nnp%n_rad(i) - 1
1581 loc = j
1582 DO k = j + 1, nnp%n_rad(i)
1583 IF (nnp%rad(i)%funccut(loc) == nnp%rad(i)%funccut(k) .AND. &
1584 nnp%rad(i)%eta(loc) == nnp%rad(i)%eta(k) .AND. &
1585 nnp%rad(i)%rs(loc) == nnp%rad(i)%rs(k) .AND. &
1586 nnp%rad(i)%nuc_ele(loc) > nnp%rad(i)%nuc_ele(k)) THEN
1587 loc = k
1588 END IF
1589 END DO
1590 CALL nnp_swaprad(nnp%rad(i), j, loc)
1591 END DO
1592
1593 DO j = 1, nnp%n_ang(i) - 1
1594 loc = j
1595 DO k = j + 1, nnp%n_ang(i)
1596 IF (nnp%ang(i)%funccut(loc) > nnp%ang(i)%funccut(k)) THEN
1597 loc = k
1598 END IF
1599 END DO
1600 CALL nnp_swapang(nnp%ang(i), j, loc)
1601 END DO
1602
1603 DO j = 1, nnp%n_ang(i) - 1
1604 loc = j
1605 DO k = j + 1, nnp%n_ang(i)
1606 IF (nnp%ang(i)%funccut(loc) == nnp%ang(i)%funccut(k) .AND. &
1607 nnp%ang(i)%eta(loc) > nnp%ang(i)%eta(k)) THEN
1608 loc = k
1609 END IF
1610 END DO
1611 CALL nnp_swapang(nnp%ang(i), j, loc)
1612 END DO
1613
1614 DO j = 1, nnp%n_ang(i) - 1
1615 loc = j
1616 DO k = j + 1, nnp%n_ang(i)
1617 IF (nnp%ang(i)%funccut(loc) == nnp%ang(i)%funccut(k) .AND. &
1618 nnp%ang(i)%eta(loc) == nnp%ang(i)%eta(k) .AND. &
1619 nnp%ang(i)%zeta(loc) > nnp%ang(i)%zeta(k)) THEN
1620 loc = k
1621 END IF
1622 END DO
1623 CALL nnp_swapang(nnp%ang(i), j, loc)
1624 END DO
1625
1626 DO j = 1, nnp%n_ang(i) - 1
1627 loc = j
1628 DO k = j + 1, nnp%n_ang(i)
1629 IF (nnp%ang(i)%funccut(loc) == nnp%ang(i)%funccut(k) .AND. &
1630 nnp%ang(i)%eta(loc) == nnp%ang(i)%eta(k) .AND. &
1631 nnp%ang(i)%zeta(loc) == nnp%ang(i)%zeta(k) .AND. &
1632 nnp%ang(i)%lam(loc) > nnp%ang(i)%lam(k)) THEN
1633 loc = k
1634 END IF
1635 END DO
1636 CALL nnp_swapang(nnp%ang(i), j, loc)
1637 END DO
1638
1639 DO j = 1, nnp%n_ang(i) - 1
1640 loc = j
1641 DO k = j + 1, nnp%n_ang(i)
1642 IF (nnp%ang(i)%funccut(loc) == nnp%ang(i)%funccut(k) .AND. &
1643 nnp%ang(i)%eta(loc) == nnp%ang(i)%eta(k) .AND. &
1644 nnp%ang(i)%zeta(loc) == nnp%ang(i)%zeta(k) .AND. &
1645 nnp%ang(i)%lam(loc) == nnp%ang(i)%lam(k) .AND. &
1646 nnp%ang(i)%nuc_ele1(loc) > nnp%ang(i)%nuc_ele1(k)) THEN
1647 loc = k
1648 END IF
1649 END DO
1650 CALL nnp_swapang(nnp%ang(i), j, loc)
1651 END DO
1652
1653 DO j = 1, nnp%n_ang(i) - 1
1654 loc = j
1655 DO k = j + 1, nnp%n_ang(i)
1656 IF (nnp%ang(i)%funccut(loc) == nnp%ang(i)%funccut(k) .AND. &
1657 nnp%ang(i)%eta(loc) == nnp%ang(i)%eta(k) .AND. &
1658 nnp%ang(i)%zeta(loc) == nnp%ang(i)%zeta(k) .AND. &
1659 nnp%ang(i)%lam(loc) == nnp%ang(i)%lam(k) .AND. &
1660 nnp%ang(i)%nuc_ele1(loc) == nnp%ang(i)%nuc_ele1(k) .AND. &
1661 nnp%ang(i)%nuc_ele2(loc) > nnp%ang(i)%nuc_ele2(k)) THEN
1662 loc = k
1663 END IF
1664 END DO
1665 CALL nnp_swapang(nnp%ang(i), j, loc)
1666 END DO
1667 END DO
1668
1669 END SUBROUTINE nnp_sort_acsf
1670
1671! **************************************************************************************************
1672!> \brief Swap two radial symmetry functions
1673!> \param rad ...
1674!> \param i ...
1675!> \param j ...
1676!> \date 2020-10-10
1677!> \author Christoph Schran (christoph.schran@rub.de)
1678! **************************************************************************************************
1679 SUBROUTINE nnp_swaprad(rad, i, j)
1680 TYPE(nnp_acsf_rad_type), INTENT(INOUT) :: rad
1681 INTEGER, INTENT(IN) :: i, j
1682
1683 CHARACTER(len=2) :: tmpc
1684 INTEGER :: tmpi
1685 REAL(kind=dp) :: tmpr
1686
1687 tmpr = rad%funccut(i)
1688 rad%funccut(i) = rad%funccut(j)
1689 rad%funccut(j) = tmpr
1690
1691 tmpr = rad%eta(i)
1692 rad%eta(i) = rad%eta(j)
1693 rad%eta(j) = tmpr
1694
1695 tmpr = rad%rs(i)
1696 rad%rs(i) = rad%rs(j)
1697 rad%rs(j) = tmpr
1698
1699 tmpc = rad%ele(i)
1700 rad%ele(i) = rad%ele(j)
1701 rad%ele(j) = tmpc
1702
1703 tmpi = rad%nuc_ele(i)
1704 rad%nuc_ele(i) = rad%nuc_ele(j)
1705 rad%nuc_ele(j) = tmpi
1706
1707 END SUBROUTINE nnp_swaprad
1708
1709! **************************************************************************************************
1710!> \brief Swap two angular symmetry functions
1711!> \param ang ...
1712!> \param i ...
1713!> \param j ...
1714!> \date 2020-10-10
1715!> \author Christoph Schran (christoph.schran@rub.de)
1716! **************************************************************************************************
1717 SUBROUTINE nnp_swapang(ang, i, j)
1718 TYPE(nnp_acsf_ang_type), INTENT(INOUT) :: ang
1719 INTEGER, INTENT(IN) :: i, j
1720
1721 CHARACTER(len=2) :: tmpc
1722 INTEGER :: tmpi
1723 REAL(kind=dp) :: tmpr
1724
1725 tmpr = ang%funccut(i)
1726 ang%funccut(i) = ang%funccut(j)
1727 ang%funccut(j) = tmpr
1728
1729 tmpr = ang%eta(i)
1730 ang%eta(i) = ang%eta(j)
1731 ang%eta(j) = tmpr
1732
1733 tmpr = ang%zeta(i)
1734 ang%zeta(i) = ang%zeta(j)
1735 ang%zeta(j) = tmpr
1736
1737 tmpr = ang%prefzeta(i)
1738 ang%prefzeta(i) = ang%prefzeta(j)
1739 ang%prefzeta(j) = tmpr
1740
1741 tmpr = ang%lam(i)
1742 ang%lam(i) = ang%lam(j)
1743 ang%lam(j) = tmpr
1744
1745 tmpc = ang%ele1(i)
1746 ang%ele1(i) = ang%ele1(j)
1747 ang%ele1(j) = tmpc
1748
1749 tmpi = ang%nuc_ele1(i)
1750 ang%nuc_ele1(i) = ang%nuc_ele1(j)
1751 ang%nuc_ele1(j) = tmpi
1752
1753 tmpc = ang%ele2(i)
1754 ang%ele2(i) = ang%ele2(j)
1755 ang%ele2(j) = tmpc
1756
1757 tmpi = ang%nuc_ele2(i)
1758 ang%nuc_ele2(i) = ang%nuc_ele2(j)
1759 ang%nuc_ele2(j) = tmpi
1760
1761 END SUBROUTINE nnp_swapang
1762
1763! **************************************************************************************************
1764!> \brief Pack symmetry functions into groups that share input parameters.
1765!> Builds nnp%rad(i)%symfgrp(:) and nnp%ang(i)%symfgrp(:) so that all
1766!> radial / angular SFs with identical (eta, rcut [, lambda, zeta])
1767!> live in one group; downstream the descriptor pass evaluates each
1768!> group's shared cutoff/exponent once and applies it across the
1769!> per-element neighbour slab.
1770!> \param nnp NNP environment with rad/ang SF metadata already populated by nnp_init_model.
1771!> \author Christoph Schran (christoph.schran@rub.de)
1772! **************************************************************************************************
1773 SUBROUTINE nnp_init_acsf_groups(nnp)
1774
1775 TYPE(nnp_type), INTENT(INOUT) :: nnp
1776
1777 INTEGER :: ang, i, izeta_tmp, j, k, m, n_symf, rad, &
1778 s, sf
1779 REAL(kind=dp) :: eta_tmp, funccut, zeta_tmp
1780
1781 DO i = 1, nnp%n_ele
1782 nnp%rad(i)%n_symfgrp = 0
1783 nnp%ang(i)%n_symfgrp = 0
1784 DO j = 1, nnp%n_ele
1785 funccut = -1.0_dp
1786 DO s = 1, nnp%n_rad(i)
1787 IF (nnp%rad(i)%ele(s) == nnp%ele(j)) THEN
1788 IF (abs(nnp%rad(i)%funccut(s) - funccut) > cutoff_eq_tol) THEN
1789 nnp%rad(i)%n_symfgrp = nnp%rad(i)%n_symfgrp + 1
1790 funccut = nnp%rad(i)%funccut(s)
1791 END IF
1792 END IF
1793 END DO
1794 END DO
1795 DO j = 1, nnp%n_ele
1796 DO k = j, nnp%n_ele
1797 funccut = -1.0_dp
1798 DO s = 1, nnp%n_ang(i)
1799 IF ((nnp%ang(i)%ele1(s) == nnp%ele(j) .AND. &
1800 nnp%ang(i)%ele2(s) == nnp%ele(k)) .OR. &
1801 (nnp%ang(i)%ele1(s) == nnp%ele(k) .AND. &
1802 nnp%ang(i)%ele2(s) == nnp%ele(j))) THEN
1803 IF (abs(nnp%ang(i)%funccut(s) - funccut) > cutoff_eq_tol) THEN
1804 nnp%ang(i)%n_symfgrp = nnp%ang(i)%n_symfgrp + 1
1805 funccut = nnp%ang(i)%funccut(s)
1806 END IF
1807 END IF
1808 END DO
1809 END DO
1810 END DO
1811 END DO
1812
1813 DO i = 1, nnp%n_ele
1814 ALLOCATE (nnp%rad(i)%symfgrp(nnp%rad(i)%n_symfgrp))
1815 ALLOCATE (nnp%ang(i)%symfgrp(nnp%ang(i)%n_symfgrp))
1816 DO j = 1, nnp%rad(i)%n_symfgrp
1817 nnp%rad(i)%symfgrp(j)%n_symf = 0
1818 END DO
1819 DO j = 1, nnp%ang(i)%n_symfgrp
1820 nnp%ang(i)%symfgrp(j)%n_symf = 0
1821 END DO
1822 END DO
1823
1824 DO i = 1, nnp%n_ele
1825 rad = 0
1826 ang = 0
1827 DO j = 1, nnp%n_ele
1828 funccut = -1.0_dp
1829 DO s = 1, nnp%n_rad(i)
1830 IF (nnp%rad(i)%ele(s) == nnp%ele(j)) THEN
1831 IF (abs(nnp%rad(i)%funccut(s) - funccut) > cutoff_eq_tol) THEN
1832 rad = rad + 1
1833 funccut = nnp%rad(i)%funccut(s)
1834 nnp%rad(i)%symfgrp(rad)%cutoff = funccut
1835 ALLOCATE (nnp%rad(i)%symfgrp(rad)%ele(1))
1836 ALLOCATE (nnp%rad(i)%symfgrp(rad)%ele_ind(1))
1837 nnp%rad(i)%symfgrp(rad)%ele(1) = nnp%ele(j)
1838 nnp%rad(i)%symfgrp(rad)%ele_ind(1) = j
1839 END IF
1840 nnp%rad(i)%symfgrp(rad)%n_symf = nnp%rad(i)%symfgrp(rad)%n_symf + 1
1841 END IF
1842 END DO
1843 END DO
1844 DO j = 1, nnp%n_ele
1845 DO k = j, nnp%n_ele
1846 funccut = -1.0_dp
1847 DO s = 1, nnp%n_ang(i)
1848 IF ((nnp%ang(i)%ele1(s) == nnp%ele(j) .AND. &
1849 nnp%ang(i)%ele2(s) == nnp%ele(k)) .OR. &
1850 (nnp%ang(i)%ele1(s) == nnp%ele(k) .AND. &
1851 nnp%ang(i)%ele2(s) == nnp%ele(j))) THEN
1852 IF (abs(nnp%ang(i)%funccut(s) - funccut) > cutoff_eq_tol) THEN
1853 ang = ang + 1
1854 funccut = nnp%ang(i)%funccut(s)
1855 nnp%ang(i)%symfgrp(ang)%cutoff = funccut
1856 ALLOCATE (nnp%ang(i)%symfgrp(ang)%ele(2))
1857 ALLOCATE (nnp%ang(i)%symfgrp(ang)%ele_ind(2))
1858 nnp%ang(i)%symfgrp(ang)%ele(1) = nnp%ele(j)
1859 nnp%ang(i)%symfgrp(ang)%ele(2) = nnp%ele(k)
1860 nnp%ang(i)%symfgrp(ang)%ele_ind(1) = j
1861 nnp%ang(i)%symfgrp(ang)%ele_ind(2) = k
1862 END IF
1863 nnp%ang(i)%symfgrp(ang)%n_symf = nnp%ang(i)%symfgrp(ang)%n_symf + 1
1864 END IF
1865 END DO
1866 END DO
1867 END DO
1868 END DO
1869
1870 DO i = 1, nnp%n_ele
1871 DO j = 1, nnp%rad(i)%n_symfgrp
1872 ALLOCATE (nnp%rad(i)%symfgrp(j)%symf(nnp%rad(i)%symfgrp(j)%n_symf))
1873 rad = 0
1874 DO s = 1, nnp%n_rad(i)
1875 IF (nnp%rad(i)%ele(s) == nnp%rad(i)%symfgrp(j)%ele(1)) THEN
1876 IF (abs(nnp%rad(i)%funccut(s) - nnp%rad(i)%symfgrp(j)%cutoff) <= cutoff_eq_tol) THEN
1877 rad = rad + 1
1878 nnp%rad(i)%symfgrp(j)%symf(rad) = s
1879 END IF
1880 END IF
1881 END DO
1882 END DO
1883 DO j = 1, nnp%ang(i)%n_symfgrp
1884 ALLOCATE (nnp%ang(i)%symfgrp(j)%symf(nnp%ang(i)%symfgrp(j)%n_symf))
1885 ang = 0
1886 DO s = 1, nnp%n_ang(i)
1887 IF ((nnp%ang(i)%ele1(s) == nnp%ang(i)%symfgrp(j)%ele(1) .AND. &
1888 nnp%ang(i)%ele2(s) == nnp%ang(i)%symfgrp(j)%ele(2)) .OR. &
1889 (nnp%ang(i)%ele1(s) == nnp%ang(i)%symfgrp(j)%ele(2) .AND. &
1890 nnp%ang(i)%ele2(s) == nnp%ang(i)%symfgrp(j)%ele(1))) THEN
1891 IF (abs(nnp%ang(i)%funccut(s) - nnp%ang(i)%symfgrp(j)%cutoff) <= cutoff_eq_tol) THEN
1892 ang = ang + 1
1893 nnp%ang(i)%symfgrp(j)%symf(ang) = s
1894 END IF
1895 END IF
1896 END DO
1897 END DO
1898 END DO
1899
1900 ! Populate packed parameter arrays for the angular SF group inner loop
1901 ! so it reads contiguous memory rather than chasing
1902 ! ang(i)%{eta,zeta,lam,prefzeta}(symf(sf)) on every iteration.
1903 ! pack_use_int_zeta and pack_izeta let the inner loop skip NINT and
1904 ! dispatch a tight integer power when zeta is integral. The radial
1905 ! group does not need packed parameters since nnp_calc_rad reads
1906 ! pretabulated splines.
1907 DO i = 1, nnp%n_ele
1908 DO j = 1, nnp%ang(i)%n_symfgrp
1909 n_symf = nnp%ang(i)%symfgrp(j)%n_symf
1910 ALLOCATE (nnp%ang(i)%symfgrp(j)%pack_eta(n_symf))
1911 ALLOCATE (nnp%ang(i)%symfgrp(j)%pack_zeta(n_symf))
1912 ALLOCATE (nnp%ang(i)%symfgrp(j)%pack_lam(n_symf))
1913 ALLOCATE (nnp%ang(i)%symfgrp(j)%pack_prefzeta(n_symf))
1914 ALLOCATE (nnp%ang(i)%symfgrp(j)%pack_izeta(n_symf))
1915 ALLOCATE (nnp%ang(i)%symfgrp(j)%pack_use_int_zeta(n_symf))
1916 DO sf = 1, n_symf
1917 m = nnp%ang(i)%symfgrp(j)%symf(sf)
1918 eta_tmp = nnp%ang(i)%eta(m)
1919 zeta_tmp = nnp%ang(i)%zeta(m)
1920 nnp%ang(i)%symfgrp(j)%pack_eta(sf) = eta_tmp
1921 nnp%ang(i)%symfgrp(j)%pack_zeta(sf) = zeta_tmp
1922 nnp%ang(i)%symfgrp(j)%pack_lam(sf) = nnp%ang(i)%lam(m)
1923 nnp%ang(i)%symfgrp(j)%pack_prefzeta(sf) = nnp%ang(i)%prefzeta(m)
1924 izeta_tmp = nint(zeta_tmp)
1925 nnp%ang(i)%symfgrp(j)%pack_izeta(sf) = izeta_tmp
1926 nnp%ang(i)%symfgrp(j)%pack_use_int_zeta(sf) = &
1927 (real(izeta_tmp, dp) == zeta_tmp)
1928 END DO
1929 END DO
1930 END DO
1931
1932 END SUBROUTINE nnp_init_acsf_groups
1933
1934! **************************************************************************************************
1935!> \brief Print a summary of the active symmetry-function set on the source rank.
1936!> Emits one line per element listing the per-element n_rad / n_ang counts
1937!> and the per-group breakdown; used at NNP init time and after sort/group
1938!> passes for traceability.
1939!> \param nnp NNP environment whose SF metadata is to be printed.
1940!> \param para_env parallel environment; only the source rank emits output.
1941!> \param printtag log-line prefix (typically "NNP" or "HELIUM_NNP").
1942!> \author Christoph Schran (christoph.schran@rub.de)
1943! **************************************************************************************************
1944 SUBROUTINE nnp_write_acsf(nnp, para_env, printtag)
1945 TYPE(nnp_type), INTENT(INOUT) :: nnp
1946 TYPE(mp_para_env_type), POINTER :: para_env
1947 CHARACTER(LEN=*), INTENT(IN) :: printtag
1948
1949 CHARACTER(len=default_string_length) :: my_label
1950 INTEGER :: i, j, unit_nr
1951 TYPE(cp_logger_type), POINTER :: logger
1952
1953 NULLIFY (logger)
1954 logger => cp_get_default_logger()
1955
1956 my_label = trim(printtag)//"| "
1957 IF (para_env%is_source()) THEN
1958 unit_nr = cp_logger_get_default_unit_nr(logger)
1959 WRITE (unit_nr, '(1X,A,1X,10(I2,1X))') trim(my_label)//" Activation functions:", nnp%actfnct(:)
1960 DO i = 1, nnp%n_ele
1961 WRITE (unit_nr, *) trim(my_label)//" short range atomic symmetry functions element "// &
1962 nnp%ele(i)//":"
1963 DO j = 1, nnp%n_rad(i)
1964 WRITE (unit_nr, '(1X,A,1X,I3,1X,A2,1X,I2,1X,A2,11X,3(F6.3,1X))') trim(my_label), j, nnp%ele(i), 2, &
1965 nnp%rad(i)%ele(j), nnp%rad(i)%eta(j), &
1966 nnp%rad(i)%rs(j), nnp%rad(i)%funccut(j)
1967 END DO
1968 DO j = 1, nnp%n_ang(i)
1969 WRITE (unit_nr, '(1X,A,1X,I3,1X,A2,1X,I2,2(1X,A2),1X,4(F6.3,1X))') &
1970 trim(my_label), j, nnp%ele(i), 3, &
1971 nnp%ang(i)%ele1(j), nnp%ang(i)%ele2(j), &
1972 nnp%ang(i)%eta(j), nnp%ang(i)%lam(j), &
1973 nnp%ang(i)%zeta(j), nnp%ang(i)%funccut(j)
1974 END DO
1975 END DO
1976 END IF
1977
1978 END SUBROUTINE nnp_write_acsf
1979
1980END MODULE nnp_acsf
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
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
Definition cp_log_handling.F:576
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
nnp_acsf
Functionality for atom centered symmetry functions for neural network potentials.
Definition nnp_acsf.F:15
nnp_acsf::nnp_calc_acsf
subroutine, public nnp_calc_acsf(nnp, i, calc_forces, stress)
Calculate atom centered symmetry functions for given atom i.
Definition nnp_acsf.F:84
nnp_acsf::nnp_write_acsf
subroutine, public nnp_write_acsf(nnp, para_env, printtag)
Print a summary of the active symmetry-function set on the source rank. Emits one line per element li...
Definition nnp_acsf.F:1945
nnp_acsf::nnp_init_acsf_groups
subroutine, public nnp_init_acsf_groups(nnp)
Pack symmetry functions into groups that share input parameters. Builds nnprad(i)symfgrp(:) and nnpan...
Definition nnp_acsf.F:1774
nnp_acsf::nnp_prepare_neighbor_cache
subroutine, public nnp_prepare_neighbor_cache(nnp)
Prepare or update the linked-cell / Verlet cache for the current geometry. Lazily allocates nnpcell_l...
Definition nnp_acsf.F:879
nnp_acsf::nnp_sort_acsf
subroutine, public nnp_sort_acsf(nnp)
Sort radial and angular symmetry functions in canonical order. Radial SFs are sorted by eta (ascendin...
Definition nnp_acsf.F:1542
nnp_acsf::nnp_sort_ele
subroutine, public nnp_sort_ele(ele, nuc_ele)
Sort an (ele, nuc_ele) pair of arrays in ascending order of atomic number. Used to canonicalise eleme...
Definition nnp_acsf.F:1503
nnp_cell_list
Linked-cell neighbour finder with a Verlet skin for the NNP descriptor. Owns the per-nnp cell-list ca...
Definition nnp_cell_list.F:18
nnp_cell_list::nnp_prepare_cell_list_cache
subroutine, public nnp_prepare_cell_list_cache(nnp)
Prepare the cell-list cache for the current force evaluation: wrap positions into the primary cell,...
Definition nnp_cell_list.F:56
nnp_cell_list::nnp_compute_neighbors_cell_list
subroutine, public nnp_compute_neighbors_cell_list(nnp, neighbor, i)
Fill the ACSF neighbour buffers for one central atom by walking the linked-cell stencil around its bi...
Definition nnp_cell_list.F:340
nnp_environment_types
Data types for neural network potentials.
Definition nnp_environment_types.F:14
nnp_environment_types::nnp_cut_tanh
integer, parameter, public nnp_cut_tanh
Definition nnp_environment_types.F:73
nnp_environment_types::nnp_cut_cos
integer, parameter, public nnp_cut_cos
derived data types
Definition nnp_environment_types.F:73
nnp_neighbor_interface
Per-nnp persistent neighbour-interface state for the NNP hot path. Separates neighbour bookkeeping fr...
Definition nnp_neighbor_interface.F:19
nnp_neighbor_interface::nnp_grp_grow_dgdr
subroutine, public nnp_grp_grow_dgdr(grp, n_needed)
Ensure a per-group dG/dr buffer holds n_needed neighbours, growing by 1.5x (with an additive floor) s...
Definition nnp_neighbor_interface.F:379
nnp_neighbor_interface::nnp_workspace_grow_caches
subroutine, public nnp_workspace_grow_caches(workspace, n_needed)
Ensure the four per-element angular cutoff caches hold at least n_needed entries. These 1D scratch ar...
Definition nnp_neighbor_interface.F:443
nnp_neighbor_interface::nnp_neighbor_interface_reset_neighbor
subroutine, public nnp_neighbor_interface_reset_neighbor(nnp, ind)
Reset per-group neighbour counters for one central element before refilling the reusable buffers....
Definition nnp_neighbor_interface.F:78
nnp_neighbor_interface::nnp_neighbor_interface_prepare
subroutine, public nnp_neighbor_interface_prepare(nnp)
Ensure pair-routing metadata and reusable workspaces are ready for the current NNP model.
Definition nnp_neighbor_interface.F:47
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::get_ptable_info
subroutine, public get_ptable_info(symbol, number, amass, ielement, covalent_radius, metallic_radius, vdw_radius, found)
Pass information about the kind given the element symbol.
Definition periodic_table.F:85
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:750
nnp_environment_types::nnp_acsf_ang_type
Set of angular symmetry function type.
Definition nnp_environment_types.F:242
nnp_environment_types::nnp_acsf_rad_type
Set of radial symmetry function type.
Definition nnp_environment_types.F:211
nnp_environment_types::nnp_symfgrp_type
Symmetry functions group type.
Definition nnp_environment_types.F:168
nnp_environment_types::nnp_type
Main data type collecting all relevant data for neural network potentials.
Definition nnp_environment_types.F:93