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