Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 : ! **************************************************************************************************
16 : MODULE 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 :
62 : CONTAINS
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 252785 : 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 252785 : CALL timeset(routineN, handle)
119 :
120 : !determine index of atom type
121 252785 : ind = nnp%ele_ind(i)
122 252785 : do_forces = calc_forces
123 :
124 : ! 1-based thread id selecting this thread's workspace column
125 252785 : tid = 1
126 252785 : !$ tid = omp_get_thread_num() + 1
127 :
128 758355 : ALLOCATE (rad_y(nnp%n_rad(ind)))
129 758355 : 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 252785 : IF (nnp%rad(ind)%n_symfgrp > 0) THEN
136 252785 : 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 252785 : n_input_nodes = nnp%neighbor_interface_state%workspace(ind, tid)%n_input_nodes
152 2487569 : 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 252785 : CALL timeset('nnp_acsf_neighbor_fill', handle_nlist)
158 1011140 : neighbor%pbc_copies = nnp%cell_list_cache%exact_pbc_copies
159 252785 : CALL nnp_neighbor_interface_reset_neighbor(nnp, ind)
160 252785 : CALL nnp_compute_neighbors_cell_list(nnp, neighbor, i)
161 252785 : CALL timestop(handle_nlist)
162 :
163 : ! Reset y:
164 252785 : rad_y = 0.0_dp
165 252785 : 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 252785 : peak = 0
171 769054 : IF (SIZE(neighbor%n_ang1) > 0) peak = MAX(peak, MAXVAL(neighbor%n_ang1))
172 769054 : IF (SIZE(neighbor%n_ang2) > 0) peak = MAX(peak, MAXVAL(neighbor%n_ang2))
173 252785 : 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 276498 : dfc_cache2 => workspace%dfc_cache2)
179 :
180 : !calc forces
181 252785 : IF (do_forces) THEN
182 : !loop over radial sym fnct grps
183 23713 : CALL timeset('nnp_acsf_radial', handle_sf)
184 74339 : DO s = 1, nnp%rad(ind)%n_symfgrp
185 50626 : n_symf_s = nnp%rad(ind)%symfgrp(s)%n_symf
186 : ! Per-group dense buffer: (3, n_symf_s, cap_s). Grown lazily.
187 50626 : CALL nnp_grp_grow_dGdr(workspace%dGdr_rad(s), neighbor%n_rad(s))
188 23713 : ASSOCIATE (rad_buf => workspace%dGdr_rad(s)%data)
189 : !loop over associated neighbors
190 1900812 : DO j = 1, neighbor%n_rad(s)
191 7400744 : rvect1 = neighbor%rad(s)%dist(1:3, j)
192 1850186 : 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 1850186 : radial_forcetmp(:, 1:n_symf_s))
196 : ! Per-group dense write: rad_buf(:, sf, j) holds dG_m/dr_j.
197 16679032 : DO sf = 1, n_symf_s
198 14778220 : m = nnp%rad(ind)%symfgrp(s)%symf(sf)
199 14778220 : self_dGdr(1, m) = self_dGdr(1, m) + radial_forcetmp(1, sf)
200 14778220 : self_dGdr(2, m) = self_dGdr(2, m) + radial_forcetmp(2, sf)
201 14778220 : self_dGdr(3, m) = self_dGdr(3, m) + radial_forcetmp(3, sf)
202 14778220 : rad_buf(1, sf, j) = -radial_forcetmp(1, sf)
203 14778220 : rad_buf(2, sf, j) = -radial_forcetmp(2, sf)
204 14778220 : rad_buf(3, sf, j) = -radial_forcetmp(3, sf)
205 14778220 : IF (PRESENT(stress)) THEN
206 7184768 : DO l = 1, 3
207 23350496 : stress(:, l, m) = stress(:, l, m) + rvect1(:)*radial_forcetmp(l, sf)
208 : END DO
209 : END IF
210 16628406 : rad_y(m) = rad_y(m) + radial_symtmp(sf)
211 : END DO
212 : END DO
213 : END ASSOCIATE
214 : END DO
215 23713 : CALL timestop(handle_sf)
216 :
217 : !loop over angular sym fnct grps
218 23713 : CALL timeset('nnp_acsf_angular', handle_sf)
219 23713 : 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 76910 : DO s = 1, nnp%ang(ind)%n_symfgrp
225 53197 : cutoff_s = nnp%ang(ind)%symfgrp(s)%cutoff
226 53197 : cutoff_sqr = cutoff_s*cutoff_s
227 53197 : n_symf_s = nnp%ang(ind)%symfgrp(s)%n_symf
228 53197 : n_ang1_s = neighbor%n_ang1(s)
229 53197 : 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 53197 : CALL nnp_grp_grow_dGdr(workspace%dGdr_ang_jj(s), n_ang1_s)
234 53197 : IF (homo_grp) THEN
235 26284 : CALL nnp_grp_grow_dGdr(workspace%dGdr_ang_kk(s), n_ang1_s)
236 26284 : n_ang2_s = 0
237 : ELSE
238 26913 : n_ang2_s = neighbor%n_ang2(s)
239 26913 : 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 51620387 : IF (n_ang1_s > 0) workspace%dGdr_ang_jj(s)%data(:, 1:n_symf_s, 1:n_ang1_s) = 0.0_dp
245 53197 : IF (homo_grp) THEN
246 20770036 : 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 15419925 : 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 53197 : 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 76910 : 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 1053340 : DO j = 1, n_ang1_s
265 4108224 : rvect1 = neighbor%ang1(s)%dist(1:3, j)
266 1027056 : r1 = neighbor%ang1(s)%dist(4, j)
267 26168785 : DO k = j + 1, n_ang1_s
268 100461780 : rvect2 = neighbor%ang1(s)%dist(1:3, k)
269 25115445 : r2 = neighbor%ang1(s)%dist(4, k)
270 25115445 : rvect3(1) = rvect2(1) - rvect1(1)
271 25115445 : rvect3(2) = rvect2(2) - rvect1(2)
272 25115445 : rvect3(3) = rvect2(3) - rvect1(3)
273 25115445 : r3_sqr = rvect3(1)*rvect3(1) + rvect3(2)*rvect3(2) + rvect3(3)*rvect3(3)
274 26142501 : IF (r3_sqr < cutoff_sqr) THEN
275 11564739 : r3 = SQRT(r3_sqr)
276 :
277 : ! -- per-triple geometry (once) --
278 11564739 : rsqr1 = r1*r1; rsqr2 = r2*r2; rsqr3 = r3*r3
279 11564739 : r2sum_il = rsqr1 + rsqr2 + rsqr3
280 11564739 : f_il = rsqr3 - rsqr1 - rsqr2
281 11564739 : g_il = -2.0_dp*r1*r2
282 11564739 : costheta_il = f_il/g_il
283 :
284 11566535 : SELECT CASE (nnp%cut_type)
285 : CASE (nnp_cut_cos)
286 1796 : arg_il = pi*r3/cutoff_s
287 1796 : fcut3_il = 0.5_dp*(COS(arg_il) + 1.0_dp)
288 1796 : dfcut3_il = -0.5_dp*SIN(arg_il)*(pi/cutoff_s)
289 : CASE (nnp_cut_tanh)
290 11562943 : tanh_il = TANH(1.0_dp - r3/cutoff_s)
291 11562943 : fcut3_il = tanh_il**3
292 11562943 : dfcut3_il = (-3.0_dp/cutoff_s)*(tanh_il**2 - tanh_il**4)
293 : CASE DEFAULT
294 11564739 : CPABORT("NNP| Cutoff function unknown")
295 : END SELECT
296 :
297 11564739 : ftot_il = fc_cache1(j)*fc_cache1(k)*fcut3_il
298 11564739 : dfcutdr1_il = dfc_cache1(j)*fc_cache1(k)*fcut3_il
299 11564739 : dfcutdr2_il = fc_cache1(j)*dfc_cache1(k)*fcut3_il
300 11564739 : dfcutdr3_il = fc_cache1(j)*fc_cache1(k)*dfcut3_il
301 :
302 11564739 : r1_inv = 1.0_dp/r1; r2_inv = 1.0_dp/r2; r3_inv = 1.0_dp/r3
303 46258956 : dr1dx_il(:) = rvect1(:)*r1_inv
304 46258956 : dr2dx_il(:) = rvect2(:)*r2_inv
305 46258956 : dr3dx_il(:) = rvect3(:)*r3_inv
306 :
307 11564739 : inv_g2_il = 1.0_dp/(g_il*g_il)
308 46258956 : DO ii = 1, 3
309 34694217 : dgdx_t1 = 2.0_dp*r2*dr1dx_il(ii)
310 34694217 : dgdx_t2 = 2.0_dp*r1*dr2dx_il(ii)
311 : dcosbase1_il(ii) = -2.0_dp*(rvect1(ii) + rvect2(ii))*g_il &
312 34694217 : - f_il*(-(dgdx_t1 + dgdx_t2))
313 : dcosbase2_il(ii) = 2.0_dp*(rvect3(ii) + rvect1(ii))*g_il &
314 34694217 : - f_il*dgdx_t1
315 : dcosbase3_il(ii) = 2.0_dp*(rvect2(ii) - rvect3(ii))*g_il &
316 46258956 : - f_il*dgdx_t2
317 : END DO
318 :
319 : ! -- fused SF loop: compute + direct scatter --
320 71451823 : DO sf = 1, n_symf_s
321 59887084 : m = off + grp_il%symf(sf)
322 59887084 : lam_il = grp_il%pack_lam(sf)
323 59887084 : zeta_il = grp_il%pack_zeta(sf)
324 59887084 : eta_il = grp_il%pack_eta(sf)
325 59887084 : prefzeta_il = grp_il%pack_prefzeta(sf)
326 :
327 59887084 : tmp_il = 1.0_dp + lam_il*costheta_il
328 59887084 : IF (tmp_il <= 0.0_dp) THEN
329 : tmpzeta_il = 0.0_dp
330 : angular_il = 0.0_dp
331 : ELSE
332 59887084 : IF (grp_il%pack_use_int_zeta(sf)) THEN
333 59887084 : izeta_il = grp_il%pack_izeta(sf)
334 59887084 : tmpzeta_il = tmp_il**(izeta_il - 1)
335 : ELSE
336 0 : tmpzeta_il = tmp_il**(zeta_il - 1.0_dp)
337 : END IF
338 59887084 : angular_il = tmpzeta_il*tmp_il
339 : END IF
340 :
341 59887084 : symtmp_il = EXP(-eta_il*r2sum_il)
342 59887084 : sym_il = prefzeta_il*angular_il*symtmp_il*ftot_il
343 59887084 : ang_y(m - off) = ang_y(m - off) + sym_il
344 :
345 59887084 : pref_lam_il = zeta_il*tmpzeta_il*lam_il*inv_g2_il
346 59887084 : tmp_il = -2.0_dp*symtmp_il*eta_il
347 59887084 : dsymdr1_il = tmp_il*r1
348 59887084 : dsymdr2_il = tmp_il*r2
349 59887084 : dsymdr3_il = tmp_il*r3
350 :
351 59887084 : pref_il = prefzeta_il*symtmp_il*ftot_il
352 59887084 : tmp1_il = prefzeta_il*angular_il*(ftot_il*dsymdr1_il + dfcutdr1_il*symtmp_il)
353 59887084 : tmp2_il = prefzeta_il*angular_il*(ftot_il*dsymdr2_il + dfcutdr2_il*symtmp_il)
354 59887084 : tmp3_il = prefzeta_il*angular_il*(ftot_il*dsymdr3_il + dfcutdr3_il*symtmp_il)
355 :
356 239548336 : DO ii = 1, 3
357 : f_jj_il(ii) = pref_il*pref_lam_il*dcosbase2_il(ii) &
358 179661252 : - tmp1_il*dr1dx_il(ii) + tmp3_il*dr3dx_il(ii)
359 : f_kk_il(ii) = pref_il*pref_lam_il*dcosbase3_il(ii) &
360 179661252 : - 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 179661252 : + tmp1_il*dr1dx_il(ii) + tmp2_il*dr2dx_il(ii)
364 179661252 : jj_buf(ii, sf, j) = jj_buf(ii, sf, j) + f_jj_il(ii)
365 239548336 : kk_buf(ii, sf, k) = kk_buf(ii, sf, k) + f_kk_il(ii)
366 : END DO
367 71451823 : IF (PRESENT(stress)) THEN
368 29217640 : DO l = 1, 3
369 : stress(:, l, m) = stress(:, l, m) &
370 94957330 : - 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 26913 : 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 1265047 : DO j = 1, n_ang1_s
388 4952536 : rvect1 = neighbor%ang1(s)%dist(1:3, j)
389 1238134 : r1 = neighbor%ang1(s)%dist(4, j)
390 44759548 : DO k = 1, n_ang2_s
391 173978004 : rvect2 = neighbor%ang2(s)%dist(1:3, k)
392 43494501 : r2 = neighbor%ang2(s)%dist(4, k)
393 43494501 : rvect3(1) = rvect2(1) - rvect1(1)
394 43494501 : rvect3(2) = rvect2(2) - rvect1(2)
395 43494501 : rvect3(3) = rvect2(3) - rvect1(3)
396 43494501 : r3_sqr = rvect3(1)*rvect3(1) + rvect3(2)*rvect3(2) + rvect3(3)*rvect3(3)
397 44732635 : IF (r3_sqr < cutoff_sqr) THEN
398 20226332 : r3 = SQRT(r3_sqr)
399 :
400 : ! -- per-triple geometry (once) --
401 20226332 : rsqr1 = r1*r1; rsqr2 = r2*r2; rsqr3 = r3*r3
402 20226332 : r2sum_il = rsqr1 + rsqr2 + rsqr3
403 20226332 : f_il = rsqr3 - rsqr1 - rsqr2
404 20226332 : g_il = -2.0_dp*r1*r2
405 20226332 : costheta_il = f_il/g_il
406 :
407 20230778 : SELECT CASE (nnp%cut_type)
408 : CASE (nnp_cut_cos)
409 4446 : arg_il = pi*r3/cutoff_s
410 4446 : fcut3_il = 0.5_dp*(COS(arg_il) + 1.0_dp)
411 4446 : dfcut3_il = -0.5_dp*SIN(arg_il)*(pi/cutoff_s)
412 : CASE (nnp_cut_tanh)
413 20221886 : tanh_il = TANH(1.0_dp - r3/cutoff_s)
414 20221886 : fcut3_il = tanh_il**3
415 20221886 : dfcut3_il = (-3.0_dp/cutoff_s)*(tanh_il**2 - tanh_il**4)
416 : CASE DEFAULT
417 20226332 : CPABORT("NNP| Cutoff function unknown")
418 : END SELECT
419 :
420 20226332 : ftot_il = fc_cache1(j)*fc_cache2(k)*fcut3_il
421 20226332 : dfcutdr1_il = dfc_cache1(j)*fc_cache2(k)*fcut3_il
422 20226332 : dfcutdr2_il = fc_cache1(j)*dfc_cache2(k)*fcut3_il
423 20226332 : dfcutdr3_il = fc_cache1(j)*fc_cache2(k)*dfcut3_il
424 :
425 20226332 : r1_inv = 1.0_dp/r1; r2_inv = 1.0_dp/r2; r3_inv = 1.0_dp/r3
426 80905328 : dr1dx_il(:) = rvect1(:)*r1_inv
427 80905328 : dr2dx_il(:) = rvect2(:)*r2_inv
428 80905328 : dr3dx_il(:) = rvect3(:)*r3_inv
429 :
430 20226332 : inv_g2_il = 1.0_dp/(g_il*g_il)
431 80905328 : DO ii = 1, 3
432 60678996 : dgdx_t1 = 2.0_dp*r2*dr1dx_il(ii)
433 60678996 : dgdx_t2 = 2.0_dp*r1*dr2dx_il(ii)
434 : dcosbase1_il(ii) = -2.0_dp*(rvect1(ii) + rvect2(ii))*g_il &
435 60678996 : - f_il*(-(dgdx_t1 + dgdx_t2))
436 : dcosbase2_il(ii) = 2.0_dp*(rvect3(ii) + rvect1(ii))*g_il &
437 60678996 : - f_il*dgdx_t1
438 : dcosbase3_il(ii) = 2.0_dp*(rvect2(ii) - rvect3(ii))*g_il &
439 80905328 : - f_il*dgdx_t2
440 : END DO
441 :
442 : ! -- fused SF loop: compute + direct scatter --
443 142017805 : DO sf = 1, n_symf_s
444 121791473 : m = off + grp_il%symf(sf)
445 121791473 : lam_il = grp_il%pack_lam(sf)
446 121791473 : zeta_il = grp_il%pack_zeta(sf)
447 121791473 : eta_il = grp_il%pack_eta(sf)
448 121791473 : prefzeta_il = grp_il%pack_prefzeta(sf)
449 :
450 121791473 : tmp_il = 1.0_dp + lam_il*costheta_il
451 121791473 : IF (tmp_il <= 0.0_dp) THEN
452 : tmpzeta_il = 0.0_dp
453 : angular_il = 0.0_dp
454 : ELSE
455 121791473 : IF (grp_il%pack_use_int_zeta(sf)) THEN
456 121791473 : izeta_il = grp_il%pack_izeta(sf)
457 121791473 : tmpzeta_il = tmp_il**(izeta_il - 1)
458 : ELSE
459 0 : tmpzeta_il = tmp_il**(zeta_il - 1.0_dp)
460 : END IF
461 121791473 : angular_il = tmpzeta_il*tmp_il
462 : END IF
463 :
464 121791473 : symtmp_il = EXP(-eta_il*r2sum_il)
465 121791473 : sym_il = prefzeta_il*angular_il*symtmp_il*ftot_il
466 121791473 : ang_y(m - off) = ang_y(m - off) + sym_il
467 :
468 121791473 : pref_lam_il = zeta_il*tmpzeta_il*lam_il*inv_g2_il
469 121791473 : tmp_il = -2.0_dp*symtmp_il*eta_il
470 121791473 : dsymdr1_il = tmp_il*r1
471 121791473 : dsymdr2_il = tmp_il*r2
472 121791473 : dsymdr3_il = tmp_il*r3
473 :
474 121791473 : pref_il = prefzeta_il*symtmp_il*ftot_il
475 121791473 : tmp1_il = prefzeta_il*angular_il*(ftot_il*dsymdr1_il + dfcutdr1_il*symtmp_il)
476 121791473 : tmp2_il = prefzeta_il*angular_il*(ftot_il*dsymdr2_il + dfcutdr2_il*symtmp_il)
477 121791473 : tmp3_il = prefzeta_il*angular_il*(ftot_il*dsymdr3_il + dfcutdr3_il*symtmp_il)
478 :
479 487165892 : DO ii = 1, 3
480 : f_jj_il(ii) = pref_il*pref_lam_il*dcosbase2_il(ii) &
481 365374419 : - tmp1_il*dr1dx_il(ii) + tmp3_il*dr3dx_il(ii)
482 : f_kk_il(ii) = pref_il*pref_lam_il*dcosbase3_il(ii) &
483 365374419 : - 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 365374419 : + tmp1_il*dr1dx_il(ii) + tmp2_il*dr2dx_il(ii)
487 365374419 : jj_buf(ii, sf, j) = jj_buf(ii, sf, j) + f_jj_il(ii)
488 487165892 : kk_buf(ii, sf, k) = kk_buf(ii, sf, k) + f_kk_il(ii)
489 : END DO
490 142017805 : IF (PRESENT(stress)) THEN
491 59415432 : DO l = 1, 3
492 : stress(:, l, m) = stress(:, l, m) &
493 193100154 : - 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 23713 : CALL timestop(handle_sf)
505 : ELSE
506 : !loop over radial sym fnct grps
507 229072 : CALL timeset('nnp_acsf_radial', handle_sf)
508 794360 : DO s = 1, nnp%rad(ind)%n_symfgrp
509 : !loop over associated neighbors
510 2724876 : DO j = 1, neighbor%n_rad(s)
511 7722064 : rvect1 = neighbor%rad(s)%dist(1:3, j)
512 1930516 : r1 = neighbor%rad(s)%dist(4, j)
513 1930516 : CALL nnp_calc_rad(nnp, ind, s, rvect1, r1, radial_symtmp(1:nnp%rad(ind)%symfgrp(s)%n_symf))
514 17187322 : DO sf = 1, nnp%rad(ind)%symfgrp(s)%n_symf
515 14691518 : m = nnp%rad(ind)%symfgrp(s)%symf(sf)
516 16622034 : rad_y(m) = rad_y(m) + radial_symtmp(sf)
517 : END DO
518 : END DO
519 : END DO
520 229072 : CALL timestop(handle_sf)
521 :
522 : !loop over angular sym fnct grps
523 229072 : CALL timeset('nnp_acsf_angular', handle_sf)
524 229072 : off = nnp%n_rad(ind)
525 692144 : DO s = 1, nnp%ang(ind)%n_symfgrp
526 463072 : cutoff_s = nnp%ang(ind)%symfgrp(s)%cutoff
527 463072 : cutoff_sqr = cutoff_s*cutoff_s
528 463072 : n_symf_s = nnp%ang(ind)%symfgrp(s)%n_symf
529 463072 : 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 463072 : nnp%cut_type, cutoff_s, fc_cache1)
534 :
535 692144 : IF (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2)) THEN
536 1121116 : DO j = 1, n_ang1_s
537 3977040 : rvect1 = neighbor%ang1(s)%dist(1:3, j)
538 994260 : r1 = neighbor%ang1(s)%dist(4, j)
539 22659329 : DO k = j + 1, n_ang1_s
540 86152852 : rvect2 = neighbor%ang1(s)%dist(1:3, k)
541 21538213 : r2 = neighbor%ang1(s)%dist(4, k)
542 21538213 : rvect3(1) = rvect2(1) - rvect1(1)
543 21538213 : rvect3(2) = rvect2(2) - rvect1(2)
544 21538213 : rvect3(3) = rvect2(3) - rvect1(3)
545 21538213 : r3_sqr = rvect3(1)*rvect3(1) + rvect3(2)*rvect3(2) + rvect3(3)*rvect3(3)
546 22532473 : IF (r3_sqr < cutoff_sqr) THEN
547 9944571 : 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 9944571 : angular_symtmp(1:n_symf_s))
551 61315369 : DO sf = 1, n_symf_s
552 51370798 : m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
553 61315369 : 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 336216 : n_ang2_s = neighbor%n_ang2(s)
561 : CALL nnp_fill_fc_cache(neighbor%ang2(s)%dist, n_ang2_s, &
562 336216 : nnp%cut_type, cutoff_s, fc_cache2)
563 :
564 1719216 : DO j = 1, n_ang1_s
565 5532000 : rvect1 = neighbor%ang1(s)%dist(1:3, j)
566 1383000 : r1 = neighbor%ang1(s)%dist(4, j)
567 39030025 : DO k = 1, n_ang2_s
568 149243236 : rvect2 = neighbor%ang2(s)%dist(1:3, k)
569 37310809 : r2 = neighbor%ang2(s)%dist(4, k)
570 37310809 : rvect3(1) = rvect2(1) - rvect1(1)
571 37310809 : rvect3(2) = rvect2(2) - rvect1(2)
572 37310809 : rvect3(3) = rvect2(3) - rvect1(3)
573 37310809 : r3_sqr = rvect3(1)*rvect3(1) + rvect3(2)*rvect3(2) + rvect3(3)*rvect3(3)
574 38693809 : IF (r3_sqr < cutoff_sqr) THEN
575 17429498 : 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 17429498 : angular_symtmp(1:n_symf_s))
579 121983157 : DO sf = 1, n_symf_s
580 104553659 : m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
581 121983157 : 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 229072 : 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 252785 : CALL nnp_check_extrapolation(nnp, ind, rad_y, ang_y)
599 :
600 252785 : IF (PRESENT(stress)) THEN
601 2112 : CALL nnp_scale_acsf(nnp, ind, do_forces, arc, rad_y, ang_y, stress)
602 : ELSE
603 250673 : CALL nnp_scale_acsf(nnp, ind, do_forces, arc, rad_y, ang_y)
604 : END IF
605 :
606 252785 : DEALLOCATE (rad_y, ang_y)
607 :
608 252785 : CALL timestop(handle)
609 :
610 252785 : 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 80110 : 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 2960552 : DO j = 1, n
632 2880442 : r_tmp = dist(4, j)
633 80110 : SELECT CASE (cut_type)
634 : CASE (nnp_cut_cos)
635 17790 : arg_tmp = pi*r_tmp/cutoff_s
636 17790 : fc_cache(j) = 0.5_dp*(COS(arg_tmp) + 1.0_dp)
637 17790 : dfc_cache(j) = -0.5_dp*SIN(arg_tmp)*(pi/cutoff_s)
638 : CASE (nnp_cut_tanh)
639 2862652 : tanh_tmp = TANH(1.0_dp - r_tmp/cutoff_s)
640 2862652 : fc_cache(j) = tanh_tmp**3
641 2880442 : dfc_cache(j) = (-3.0_dp/cutoff_s)*(tanh_tmp**2 - tanh_tmp**4)
642 : END SELECT
643 : END DO
644 :
645 80110 : 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 799288 : 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 3831208 : DO j = 1, n
666 3031920 : r_tmp = dist(4, j)
667 799288 : SELECT CASE (cut_type)
668 : CASE (nnp_cut_cos)
669 584048 : fc_cache(j) = 0.5_dp*(COS(pi*r_tmp/cutoff_s) + 1.0_dp)
670 : CASE (nnp_cut_tanh)
671 2447872 : tanh_tmp = TANH(1.0_dp - r_tmp/cutoff_s)
672 3031920 : fc_cache(j) = tanh_tmp**3
673 : END SELECT
674 : END DO
675 :
676 799288 : 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 55492 : SUBROUTINE nnp_prepare_neighbor_cache(nnp)
689 :
690 : TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
691 :
692 : INTEGER :: ind
693 :
694 56546 : IF (.NOT. ALLOCATED(nnp%cell_list_cache)) ALLOCATE (nnp%cell_list_cache)
695 55492 : IF (.NOT. ALLOCATED(nnp%neighbor_interface_state)) ALLOCATE (nnp%neighbor_interface_state)
696 55492 : CALL nnp_prepare_cell_list_cache(nnp)
697 55492 : 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 221613 : DO ind = 1, nnp%n_ele
703 221613 : IF (nnp%rad(ind)%n_symfgrp > 0) THEN
704 166138 : IF (.NOT. nnp%rad(ind)%symfgrp(1)%spline_built) THEN
705 17 : CALL nnp_build_radial_splines(nnp)
706 17 : EXIT
707 : END IF
708 : END IF
709 : END DO
710 :
711 55492 : 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 252785 : 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 252785 : extrapolate = .FALSE.
734 :
735 3800797 : DO j = 1, nnp%n_rad(ind)
736 3800797 : IF (rad_y(j) - nnp%rad(ind)%loc_max(j) > threshold) THEN
737 : extrapolate = .TRUE.
738 3548012 : ELSE IF (-rad_y(j) + nnp%rad(ind)%loc_min(j) > threshold) THEN
739 164 : extrapolate = .TRUE.
740 : END IF
741 : END DO
742 2127593 : DO j = 1, nnp%n_ang(ind)
743 2127593 : IF (ang_y(j) - nnp%ang(ind)%loc_max(j) > threshold) THEN
744 : extrapolate = .TRUE.
745 1874808 : ELSE IF (-ang_y(j) + nnp%ang(ind)%loc_min(j) > threshold) THEN
746 30 : 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 252785 : IF (extrapolate) nnp%output_expol = .TRUE.
753 :
754 252785 : 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 252785 : 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 252785 : tid = 1
784 252785 : !$ tid = omp_get_thread_num() + 1
785 :
786 : ! INOUT (not OUT): stress is per-input-node and accumulates across central atoms.
787 :
788 252785 : IF (nnp%center_acsf) THEN
789 3800797 : DO j = 1, nnp%n_rad(ind)
790 3800797 : arc%layer(1)%node(j) = rad_y(j) - nnp%rad(ind)%loc_av(j)
791 : END DO
792 252785 : off = nnp%n_rad(ind)
793 2127593 : DO j = 1, nnp%n_ang(ind)
794 2127593 : arc%layer(1)%node(j + off) = ang_y(j) - nnp%ang(ind)%loc_av(j)
795 : END DO
796 :
797 252785 : IF (nnp%scale_acsf) THEN
798 3800797 : DO j = 1, nnp%n_rad(ind)
799 : arc%layer(1)%node(j) = arc%layer(1)%node(j)/ &
800 3800797 : (nnp%rad(ind)%loc_max(j) - nnp%rad(ind)%loc_min(j))*(nnp%scmax - nnp%scmin) + nnp%scmin
801 : END DO
802 252785 : off = nnp%n_rad(ind)
803 2127593 : DO j = 1, nnp%n_ang(ind)
804 : arc%layer(1)%node(j + off) = arc%layer(1)%node(j + off)/ &
805 2127593 : (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 0 : ELSE IF (nnp%scale_acsf) THEN
809 0 : 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 0 : (nnp%scmax - nnp%scmin) + nnp%scmin
813 : END DO
814 0 : off = nnp%n_rad(ind)
815 0 : 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 0 : (nnp%scmax - nnp%scmin) + nnp%scmin
819 : END DO
820 0 : ELSE IF (nnp%scale_sigma_acsf) THEN
821 0 : DO j = 1, nnp%n_rad(ind)
822 : arc%layer(1)%node(j) = (rad_y(j) - nnp%rad(ind)%loc_av(j))/ &
823 0 : nnp%rad(ind)%sigma(j)*(nnp%scmax - nnp%scmin) + nnp%scmin
824 : END DO
825 0 : off = nnp%n_rad(ind)
826 0 : DO j = 1, nnp%n_ang(ind)
827 : arc%layer(1)%node(j + off) = (ang_y(j) - nnp%ang(ind)%loc_av(j))/ &
828 0 : nnp%ang(ind)%sigma(j)*(nnp%scmax - nnp%scmin) + nnp%scmin
829 : END DO
830 : ELSE
831 0 : DO j = 1, nnp%n_rad(ind)
832 0 : arc%layer(1)%node(j) = rad_y(j)
833 : END DO
834 0 : off = nnp%n_rad(ind)
835 0 : DO j = 1, nnp%n_ang(ind)
836 0 : arc%layer(1)%node(j + off) = ang_y(j)
837 : END DO
838 : END IF
839 :
840 252785 : 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 74339 : DO s = 1, nnp%rad(ind)%n_symfgrp
852 50626 : n_symf_s = nnp%rad(ind)%symfgrp(s)%n_symf
853 23713 : ASSOCIATE (rad_buf => workspace%dGdr_rad(s)%data)
854 415634 : DO sf = 1, n_symf_s
855 365008 : m = nnp%rad(ind)%symfgrp(s)%symf(sf)
856 365008 : IF (nnp%scale_acsf) THEN
857 : scale = (nnp%scmax - nnp%scmin)/ &
858 365008 : (nnp%rad(ind)%loc_max(m) - nnp%rad(ind)%loc_min(m))
859 : ELSE
860 0 : scale = (nnp%scmax - nnp%scmin)/nnp%rad(ind)%sigma(m)
861 : END IF
862 365008 : self_dGdr(1, m) = self_dGdr(1, m)*scale
863 365008 : self_dGdr(2, m) = self_dGdr(2, m)*scale
864 365008 : self_dGdr(3, m) = self_dGdr(3, m)*scale
865 15193854 : DO j = 1, neighbor%n_rad(s)
866 14778220 : rad_buf(1, sf, j) = rad_buf(1, sf, j)*scale
867 14778220 : rad_buf(2, sf, j) = rad_buf(2, sf, j)*scale
868 15143228 : 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 23713 : off = nnp%n_rad(ind)
876 100623 : DO s = 1, nnp%ang(ind)%n_symfgrp
877 53197 : n_symf_s = nnp%ang(ind)%symfgrp(s)%n_symf
878 53197 : n_ang1_s = neighbor%n_ang1(s)
879 53197 : homo_grp = (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2))
880 53197 : 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 26913 : n_ang2_s = neighbor%n_ang2(s)
885 : END IF
886 : ASSOCIATE (jj_buf => workspace%dGdr_ang_jj(s)%data, &
887 23713 : kk_buf => workspace%dGdr_ang_kk(s)%data)
888 304153 : DO sf = 1, n_symf_s
889 250956 : m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
890 250956 : IF (nnp%scale_acsf) THEN
891 : scale = (nnp%scmax - nnp%scmin)/ &
892 250956 : (nnp%ang(ind)%loc_max(m - off) - nnp%ang(ind)%loc_min(m - off))
893 : ELSE
894 0 : scale = (nnp%scmax - nnp%scmin)/nnp%ang(ind)%sigma(m - off)
895 : END IF
896 250956 : self_dGdr(1, m) = self_dGdr(1, m)*scale
897 250956 : self_dGdr(2, m) = self_dGdr(2, m)*scale
898 250956 : self_dGdr(3, m) = self_dGdr(3, m)*scale
899 12576456 : DO j = 1, n_ang1_s
900 12325500 : jj_buf(1, sf, j) = jj_buf(1, sf, j)*scale
901 12325500 : jj_buf(2, sf, j) = jj_buf(2, sf, j)*scale
902 12576456 : jj_buf(3, sf, j) = jj_buf(3, sf, j)*scale
903 : END DO
904 8927767 : DO k = 1, n_ang2_s
905 8623614 : kk_buf(1, sf, k) = kk_buf(1, sf, k)*scale
906 8623614 : kk_buf(2, sf, k) = kk_buf(2, sf, k)*scale
907 8874570 : 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 252785 : IF (PRESENT(stress)) THEN
917 2112 : IF (nnp%scale_acsf) THEN
918 35904 : 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 441408 : (nnp%scmax - nnp%scmin)
921 : END DO
922 2112 : off = nnp%n_rad(ind)
923 27456 : 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 331584 : (nnp%scmax - nnp%scmin)
927 : END DO
928 0 : ELSE IF (nnp%scale_sigma_acsf) THEN
929 0 : DO j = 1, nnp%n_rad(ind)
930 0 : stress(:, :, j) = stress(:, :, j)/nnp%rad(ind)%sigma(j)*(nnp%scmax - nnp%scmin)
931 : END DO
932 0 : off = nnp%n_rad(ind)
933 0 : DO j = 1, nnp%n_ang(ind)
934 0 : 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 252785 : 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 3780702 : 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 3780702 : grp => nnp%rad(ind)%symfgrp(s)
971 3780702 : 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 3780702 : IF (r >= grp%spline_x_max) THEN
982 0 : DO sf = 1, n_symf
983 0 : sym(sf) = 0.0_dp
984 : END DO
985 0 : IF (PRESENT(force)) THEN
986 0 : DO sf = 1, n_symf
987 0 : force(1, sf) = 0.0_dp
988 0 : force(2, sf) = 0.0_dp
989 0 : force(3, sf) = 0.0_dp
990 : END DO
991 : END IF
992 0 : RETURN
993 : END IF
994 :
995 3780702 : i = INT(r*grp%spline_dx_inv) + 1
996 3780702 : IF (i < 1) i = 1
997 3780702 : IF (i > grp%spline_n - 1) i = grp%spline_n - 1
998 :
999 3780702 : t = (r - REAL(i - 1, KIND=dp)*grp%spline_dx)*grp%spline_dx_inv
1000 3780702 : t2 = t*t
1001 3780702 : t3 = t2*t
1002 :
1003 3780702 : h00 = 2.0_dp*t3 - 3.0_dp*t2 + 1.0_dp
1004 3780702 : h10 = t3 - 2.0_dp*t2 + t
1005 3780702 : h01 = -2.0_dp*t3 + 3.0_dp*t2
1006 3780702 : h11 = t3 - t2
1007 3780702 : h10_dx = h10*grp%spline_dx
1008 3780702 : h11_dx = h11*grp%spline_dx
1009 :
1010 3780702 : IF (PRESENT(force)) THEN
1011 1850186 : dh00 = 6.0_dp*(t2 - t)
1012 1850186 : dh10 = 3.0_dp*t2 - 4.0_dp*t + 1.0_dp
1013 1850186 : dh01 = -dh00
1014 1850186 : dh11 = 3.0_dp*t2 - 2.0_dp*t
1015 :
1016 1850186 : r_inv = 1.0_dp/r
1017 1850186 : drdx_x = rvect(1)*r_inv
1018 1850186 : drdx_y = rvect(2)*r_inv
1019 1850186 : 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 1850186 : !$OMP SIMD PRIVATE(yi, yi1, dyi, dyi1, dsymdr_full)
1026 : DO sf = 1, n_symf
1027 14778220 : yi = spy_i(sf)
1028 14778220 : yi1 = spy_i1(sf)
1029 14778220 : dyi = spdy_i(sf)
1030 14778220 : dyi1 = spdy_i1(sf)
1031 14778220 : sym(sf) = h00*yi + h10_dx*dyi + h01*yi1 + h11_dx*dyi1
1032 14778220 : dsymdr_full = (dh00*yi + dh01*yi1)*grp%spline_dx_inv + dh10*dyi + dh11*dyi1
1033 14778220 : force(1, sf) = dsymdr_full*drdx_x
1034 14778220 : force(2, sf) = dsymdr_full*drdx_y
1035 14778220 : 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 1930516 : !$OMP SIMD
1042 : DO sf = 1, n_symf
1043 : sym(sf) = h00*spy_i(sf) + h10_dx*spdy_i(sf) + &
1044 14691518 : 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 17 : 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 17 : CALL timeset(routineN, handle)
1068 :
1069 52 : DO ind = 1, nnp%n_ele
1070 123 : DO s = 1, nnp%rad(ind)%n_symfgrp
1071 35 : ASSOCIATE (grp => nnp%rad(ind)%symfgrp(s))
1072 71 : n_symf = grp%n_symf
1073 71 : cutoff = grp%cutoff
1074 71 : dr = cutoff/REAL(nnp%rad_spline_n - 1, KIND=dp)
1075 :
1076 0 : IF (ALLOCATED(grp%spline_y)) DEALLOCATE (grp%spline_y)
1077 71 : IF (ALLOCATED(grp%spline_dy)) DEALLOCATE (grp%spline_dy)
1078 284 : ALLOCATE (grp%spline_y(MAX(1, n_symf), nnp%rad_spline_n))
1079 284 : ALLOCATE (grp%spline_dy(MAX(1, n_symf), nnp%rad_spline_n))
1080 71 : grp%spline_n = nnp%rad_spline_n
1081 71 : grp%spline_dx = dr
1082 71 : grp%spline_dx_inv = 1.0_dp/dr
1083 71 : 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 621 : DO sf = 1, n_symf
1088 550 : k = grp%symf(sf)
1089 550 : eta = nnp%rad(ind)%eta(k)
1090 550 : rs = nnp%rad(ind)%rs(k)
1091 :
1092 4506221 : DO p = 1, nnp%rad_spline_n
1093 4505600 : r = REAL(p - 1, KIND=dp)*dr
1094 :
1095 4816896 : SELECT CASE (nnp%cut_type)
1096 : CASE (nnp_cut_cos)
1097 311296 : arg = pi*r/cutoff
1098 311296 : fcut = 0.5_dp*(COS(arg) + 1.0_dp)
1099 311296 : dfcutdr = -0.5_dp*SIN(arg)*(pi/cutoff)
1100 : CASE (nnp_cut_tanh)
1101 4194304 : tanh_tmp = TANH(1.0_dp - r/cutoff)
1102 4194304 : fcut = tanh_tmp**3
1103 4194304 : dfcutdr = (-3.0_dp/cutoff)*(tanh_tmp**2 - tanh_tmp**4)
1104 : CASE DEFAULT
1105 4505600 : CPABORT("NNP| Cutoff function unknown")
1106 : END SELECT
1107 :
1108 4505600 : exp_term = EXP(-eta*(r - rs)**2)
1109 :
1110 4505600 : grp%spline_y(sf, p) = exp_term*fcut
1111 : grp%spline_dy(sf, p) = exp_term*(-2.0_dp*eta*(r - rs))*fcut + &
1112 4506150 : 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 621 : DO sf = 1, n_symf
1119 550 : grp%spline_y(sf, nnp%rad_spline_n) = 0.0_dp
1120 621 : grp%spline_dy(sf, nnp%rad_spline_n) = 0.0_dp
1121 : END DO
1122 :
1123 142 : grp%spline_built = .TRUE.
1124 : END ASSOCIATE
1125 : END DO
1126 : END DO
1127 :
1128 17 : CALL timestop(handle)
1129 :
1130 17 : 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 27374069 : SUBROUTINE nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, r1, r2, r3, &
1174 27374069 : 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 54748138 : 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 27374069 : do_forces = PRESENT(force)
1202 27374069 : grp => nnp%ang(ind)%symfgrp(s)
1203 27374069 : n_symf = grp%n_symf
1204 27374069 : fcut_rc = grp%cutoff
1205 :
1206 27374069 : rsqr1 = r1*r1
1207 27374069 : rsqr2 = r2*r2
1208 27374069 : rsqr3 = r3*r3
1209 27374069 : r2sum = rsqr1 + rsqr2 + rsqr3
1210 :
1211 27374069 : f = rsqr3 - rsqr1 - rsqr2
1212 27374069 : g = -2.0_dp*r1*r2
1213 27374069 : costheta = f/g
1214 :
1215 : ! Compute fcut3 for r3 (j-k distance -- cannot be precomputed by caller)
1216 27573631 : SELECT CASE (nnp%cut_type)
1217 : CASE (nnp_cut_cos)
1218 199562 : arg_tmp = pi*r3/fcut_rc
1219 199562 : fcut3 = 0.5_dp*(COS(arg_tmp) + 1.0_dp)
1220 199562 : IF (do_forces) dfcut3 = -0.5_dp*SIN(arg_tmp)*(pi/fcut_rc)
1221 : CASE (nnp_cut_tanh)
1222 27174507 : tanh_tmp = TANH(1.0_dp - r3/fcut_rc)
1223 27174507 : fcut3 = tanh_tmp**3
1224 27174507 : IF (do_forces) dfcut3 = (-3.0_dp/fcut_rc)*(tanh_tmp**2 - tanh_tmp**4)
1225 : CASE DEFAULT
1226 27374069 : CPABORT("NNP| Cutoff function unknown")
1227 : END SELECT
1228 :
1229 : ! Use precomputed fcut values for j and k neighbors
1230 27374069 : ftot = fcut_j*fcut_k*fcut3
1231 :
1232 27374069 : IF (do_forces) THEN
1233 : ! Combined cutoff derivatives (product rule)
1234 0 : dfcutdr1 = dfcut_j*fcut_k*fcut3
1235 0 : dfcutdr2 = fcut_j*dfcut_k*fcut3
1236 0 : dfcutdr3 = fcut_j*fcut_k*dfcut3
1237 :
1238 0 : dr1dx(:) = rvect1(:)/r1
1239 0 : dr2dx(:) = rvect2(:)/r2
1240 0 : 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 0 : inv_g2 = 1.0_dp/(g*g)
1247 0 : DO ii = 1, 3
1248 0 : tmp1 = 2.0_dp*r2*dr1dx(ii)
1249 0 : tmp2 = 2.0_dp*r1*dr2dx(ii)
1250 0 : dgdx1(ii) = -(tmp1 + tmp2)
1251 : dgdx2(ii) = tmp1
1252 : dgdx3(ii) = tmp2
1253 :
1254 0 : dcosbase1(ii) = -2.0_dp*(rvect1(ii) + rvect2(ii))*g - f*dgdx1(ii)
1255 0 : dcosbase2(ii) = 2.0_dp*(rvect3(ii) + rvect1(ii))*g - f*dgdx2(ii)
1256 0 : 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 183298526 : DO sf = 1, n_symf
1264 155924457 : tmp = 1.0_dp + grp%pack_lam(sf)*costheta
1265 183298526 : IF (tmp <= 0.0_dp) THEN
1266 0 : tmpzeta_arr(sf) = 0.0_dp
1267 0 : angular_arr(sf) = 0.0_dp
1268 : ELSE
1269 155924457 : IF (grp%pack_use_int_zeta(sf)) THEN
1270 155924457 : izeta = grp%pack_izeta(sf)
1271 155924457 : tmpzeta_arr(sf) = tmp**(izeta - 1)
1272 : ELSE
1273 0 : tmpzeta_arr(sf) = tmp**(grp%pack_zeta(sf) - 1.0_dp)
1274 : END IF
1275 155924457 : 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 183298526 : !$OMP SIMD
1286 : DO sf = 1, n_symf
1287 155924457 : symtmp_arr(sf) = EXP(-grp%pack_eta(sf)*r2sum)
1288 : END DO
1289 :
1290 : ! ---- Pass 3: sym scatter (vectorizable) ----
1291 27374069 : !$OMP SIMD
1292 : DO sf = 1, n_symf
1293 155924457 : 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 27374069 : IF (do_forces) THEN
1298 0 : DO sf = 1, n_symf
1299 0 : symtmp = symtmp_arr(sf)
1300 0 : angular = angular_arr(sf)
1301 0 : tmpzeta = tmpzeta_arr(sf)
1302 0 : eta = grp%pack_eta(sf)
1303 0 : lam = grp%pack_lam(sf)
1304 0 : zeta = grp%pack_zeta(sf)
1305 0 : 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 0 : pref_lam = zeta*tmpzeta*lam*inv_g2
1310 :
1311 0 : tmp = -2.0_dp*symtmp*eta
1312 0 : dsymdr1 = tmp*r1
1313 0 : dsymdr2 = tmp*r2
1314 0 : dsymdr3 = tmp*r3
1315 :
1316 0 : tmp = prefzeta*symtmp*ftot
1317 0 : tmp1 = prefzeta*angular*(ftot*dsymdr1 + dfcutdr1*symtmp)
1318 0 : tmp2 = prefzeta*angular*(ftot*dsymdr2 + dfcutdr2*symtmp)
1319 0 : tmp3 = prefzeta*angular*(ftot*dsymdr3 + dfcutdr3*symtmp)
1320 0 : DO ii = 1, 3
1321 0 : force(ii, 1, sf) = tmp*pref_lam*dcosbase1(ii) + tmp1*dr1dx(ii) + tmp2*dr2dx(ii)
1322 0 : force(ii, 2, sf) = tmp*pref_lam*dcosbase2(ii) - tmp1*dr1dx(ii) + tmp3*dr3dx(ii)
1323 0 : 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 27374069 : 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 17 : 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 52 : DO i = 1, SIZE(ele)
1347 52 : CALL get_ptable_info(ele(i), number=nuc_ele(i))
1348 : END DO
1349 :
1350 35 : DO i = 1, SIZE(ele) - 1
1351 18 : minimum = nuc_ele(i)
1352 18 : loc = i
1353 37 : DO j = i + 1, SIZE(ele)
1354 37 : IF (nuc_ele(j) < minimum) THEN
1355 18 : loc = j
1356 18 : minimum = nuc_ele(j)
1357 : END IF
1358 : END DO
1359 18 : tmp_nuc_ele = nuc_ele(i)
1360 18 : nuc_ele(i) = nuc_ele(loc)
1361 18 : nuc_ele(loc) = tmp_nuc_ele
1362 :
1363 18 : tmp_ele = ele(i)
1364 18 : ele(i) = ele(loc)
1365 35 : ele(loc) = tmp_ele
1366 : END DO
1367 :
1368 17 : 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 17 : SUBROUTINE nnp_sort_acsf(nnp)
1379 : TYPE(nnp_type), INTENT(INOUT) :: nnp
1380 :
1381 : INTEGER :: i, j, k, loc
1382 :
1383 52 : DO i = 1, nnp%n_ele
1384 550 : DO j = 1, nnp%n_rad(i) - 1
1385 515 : loc = j
1386 4591 : DO k = j + 1, nnp%n_rad(i)
1387 4591 : IF (nnp%rad(i)%funccut(loc) > nnp%rad(i)%funccut(k)) THEN
1388 6 : loc = k
1389 : END IF
1390 : END DO
1391 550 : CALL nnp_swaprad(nnp%rad(i), j, loc)
1392 : END DO
1393 :
1394 550 : DO j = 1, nnp%n_rad(i) - 1
1395 515 : loc = j
1396 4591 : DO k = j + 1, nnp%n_rad(i)
1397 4076 : IF (nnp%rad(i)%funccut(loc) == nnp%rad(i)%funccut(k) .AND. &
1398 515 : nnp%rad(i)%eta(loc) > nnp%rad(i)%eta(k)) THEN
1399 552 : loc = k
1400 : END IF
1401 : END DO
1402 550 : CALL nnp_swaprad(nnp%rad(i), j, loc)
1403 : END DO
1404 :
1405 550 : DO j = 1, nnp%n_rad(i) - 1
1406 515 : loc = j
1407 4591 : DO k = j + 1, nnp%n_rad(i)
1408 : IF (nnp%rad(i)%funccut(loc) == nnp%rad(i)%funccut(k) .AND. &
1409 4076 : nnp%rad(i)%eta(loc) == nnp%rad(i)%eta(k) .AND. &
1410 515 : nnp%rad(i)%rs(loc) > nnp%rad(i)%rs(k)) THEN
1411 64 : loc = k
1412 : END IF
1413 : END DO
1414 550 : CALL nnp_swaprad(nnp%rad(i), j, loc)
1415 : END DO
1416 :
1417 550 : DO j = 1, nnp%n_rad(i) - 1
1418 515 : loc = j
1419 4591 : 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 4076 : nnp%rad(i)%rs(loc) == nnp%rad(i)%rs(k) .AND. &
1423 515 : nnp%rad(i)%nuc_ele(loc) > nnp%rad(i)%nuc_ele(k)) THEN
1424 0 : loc = k
1425 : END IF
1426 : END DO
1427 550 : CALL nnp_swaprad(nnp%rad(i), j, loc)
1428 : END DO
1429 :
1430 420 : DO j = 1, nnp%n_ang(i) - 1
1431 385 : loc = j
1432 2778 : DO k = j + 1, nnp%n_ang(i)
1433 2778 : IF (nnp%ang(i)%funccut(loc) > nnp%ang(i)%funccut(k)) THEN
1434 3 : loc = k
1435 : END IF
1436 : END DO
1437 420 : CALL nnp_swapang(nnp%ang(i), j, loc)
1438 : END DO
1439 :
1440 420 : DO j = 1, nnp%n_ang(i) - 1
1441 385 : loc = j
1442 2778 : DO k = j + 1, nnp%n_ang(i)
1443 2393 : IF (nnp%ang(i)%funccut(loc) == nnp%ang(i)%funccut(k) .AND. &
1444 385 : nnp%ang(i)%eta(loc) > nnp%ang(i)%eta(k)) THEN
1445 448 : loc = k
1446 : END IF
1447 : END DO
1448 420 : CALL nnp_swapang(nnp%ang(i), j, loc)
1449 : END DO
1450 :
1451 420 : DO j = 1, nnp%n_ang(i) - 1
1452 385 : loc = j
1453 2778 : DO k = j + 1, nnp%n_ang(i)
1454 : IF (nnp%ang(i)%funccut(loc) == nnp%ang(i)%funccut(k) .AND. &
1455 2393 : nnp%ang(i)%eta(loc) == nnp%ang(i)%eta(k) .AND. &
1456 385 : nnp%ang(i)%zeta(loc) > nnp%ang(i)%zeta(k)) THEN
1457 7 : loc = k
1458 : END IF
1459 : END DO
1460 420 : CALL nnp_swapang(nnp%ang(i), j, loc)
1461 : END DO
1462 :
1463 420 : DO j = 1, nnp%n_ang(i) - 1
1464 385 : loc = j
1465 2778 : 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 2393 : nnp%ang(i)%zeta(loc) == nnp%ang(i)%zeta(k) .AND. &
1469 385 : nnp%ang(i)%lam(loc) > nnp%ang(i)%lam(k)) THEN
1470 168 : loc = k
1471 : END IF
1472 : END DO
1473 420 : CALL nnp_swapang(nnp%ang(i), j, loc)
1474 : END DO
1475 :
1476 420 : DO j = 1, nnp%n_ang(i) - 1
1477 385 : loc = j
1478 2778 : 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 2393 : nnp%ang(i)%lam(loc) == nnp%ang(i)%lam(k) .AND. &
1483 385 : nnp%ang(i)%nuc_ele1(loc) > nnp%ang(i)%nuc_ele1(k)) THEN
1484 48 : loc = k
1485 : END IF
1486 : END DO
1487 420 : CALL nnp_swapang(nnp%ang(i), j, loc)
1488 : END DO
1489 :
1490 437 : DO j = 1, nnp%n_ang(i) - 1
1491 385 : loc = j
1492 2778 : 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 2393 : nnp%ang(i)%nuc_ele1(loc) == nnp%ang(i)%nuc_ele1(k) .AND. &
1498 385 : nnp%ang(i)%nuc_ele2(loc) > nnp%ang(i)%nuc_ele2(k)) THEN
1499 33 : loc = k
1500 : END IF
1501 : END DO
1502 420 : CALL nnp_swapang(nnp%ang(i), j, loc)
1503 : END DO
1504 : END DO
1505 :
1506 17 : 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 2060 : 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 2060 : tmpr = rad%funccut(i)
1525 2060 : rad%funccut(i) = rad%funccut(j)
1526 2060 : rad%funccut(j) = tmpr
1527 :
1528 2060 : tmpr = rad%eta(i)
1529 2060 : rad%eta(i) = rad%eta(j)
1530 2060 : rad%eta(j) = tmpr
1531 :
1532 2060 : tmpr = rad%rs(i)
1533 2060 : rad%rs(i) = rad%rs(j)
1534 2060 : rad%rs(j) = tmpr
1535 :
1536 2060 : tmpc = rad%ele(i)
1537 2060 : rad%ele(i) = rad%ele(j)
1538 2060 : rad%ele(j) = tmpc
1539 :
1540 2060 : tmpi = rad%nuc_ele(i)
1541 2060 : rad%nuc_ele(i) = rad%nuc_ele(j)
1542 2060 : rad%nuc_ele(j) = tmpi
1543 :
1544 2060 : 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 2310 : 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 2310 : tmpr = ang%funccut(i)
1563 2310 : ang%funccut(i) = ang%funccut(j)
1564 2310 : ang%funccut(j) = tmpr
1565 :
1566 2310 : tmpr = ang%eta(i)
1567 2310 : ang%eta(i) = ang%eta(j)
1568 2310 : ang%eta(j) = tmpr
1569 :
1570 2310 : tmpr = ang%zeta(i)
1571 2310 : ang%zeta(i) = ang%zeta(j)
1572 2310 : ang%zeta(j) = tmpr
1573 :
1574 2310 : tmpr = ang%prefzeta(i)
1575 2310 : ang%prefzeta(i) = ang%prefzeta(j)
1576 2310 : ang%prefzeta(j) = tmpr
1577 :
1578 2310 : tmpr = ang%lam(i)
1579 2310 : ang%lam(i) = ang%lam(j)
1580 2310 : ang%lam(j) = tmpr
1581 :
1582 2310 : tmpc = ang%ele1(i)
1583 2310 : ang%ele1(i) = ang%ele1(j)
1584 2310 : ang%ele1(j) = tmpc
1585 :
1586 2310 : tmpi = ang%nuc_ele1(i)
1587 2310 : ang%nuc_ele1(i) = ang%nuc_ele1(j)
1588 2310 : ang%nuc_ele1(j) = tmpi
1589 :
1590 2310 : tmpc = ang%ele2(i)
1591 2310 : ang%ele2(i) = ang%ele2(j)
1592 2310 : ang%ele2(j) = tmpc
1593 :
1594 2310 : tmpi = ang%nuc_ele2(i)
1595 2310 : ang%nuc_ele2(i) = ang%nuc_ele2(j)
1596 2310 : ang%nuc_ele2(j) = tmpi
1597 :
1598 2310 : 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 17 : 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 52 : DO i = 1, nnp%n_ele
1619 35 : nnp%rad(i)%n_symfgrp = 0
1620 35 : nnp%ang(i)%n_symfgrp = 0
1621 108 : DO j = 1, nnp%n_ele
1622 73 : funccut = -1.0_dp
1623 1246 : DO s = 1, nnp%n_rad(i)
1624 1211 : IF (nnp%rad(i)%ele(s) == nnp%ele(j)) THEN
1625 550 : IF (ABS(nnp%rad(i)%funccut(s) - funccut) > cutoff_eq_tol) THEN
1626 71 : nnp%rad(i)%n_symfgrp = nnp%rad(i)%n_symfgrp + 1
1627 71 : funccut = nnp%rad(i)%funccut(s)
1628 : END IF
1629 : END IF
1630 : END DO
1631 : END DO
1632 125 : DO j = 1, nnp%n_ele
1633 222 : DO k = j, nnp%n_ele
1634 114 : funccut = -1.0_dp
1635 1507 : DO s = 1, nnp%n_ang(i)
1636 : IF ((nnp%ang(i)%ele1(s) == nnp%ele(j) .AND. &
1637 1320 : nnp%ang(i)%ele2(s) == nnp%ele(k)) .OR. &
1638 : (nnp%ang(i)%ele1(s) == nnp%ele(k) .AND. &
1639 114 : nnp%ang(i)%ele2(s) == nnp%ele(j))) THEN
1640 420 : IF (ABS(nnp%ang(i)%funccut(s) - funccut) > cutoff_eq_tol) THEN
1641 86 : nnp%ang(i)%n_symfgrp = nnp%ang(i)%n_symfgrp + 1
1642 86 : 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 52 : DO i = 1, nnp%n_ele
1651 176 : ALLOCATE (nnp%rad(i)%symfgrp(nnp%rad(i)%n_symfgrp))
1652 191 : ALLOCATE (nnp%ang(i)%symfgrp(nnp%ang(i)%n_symfgrp))
1653 106 : DO j = 1, nnp%rad(i)%n_symfgrp
1654 106 : nnp%rad(i)%symfgrp(j)%n_symf = 0
1655 : END DO
1656 138 : DO j = 1, nnp%ang(i)%n_symfgrp
1657 121 : nnp%ang(i)%symfgrp(j)%n_symf = 0
1658 : END DO
1659 : END DO
1660 :
1661 52 : DO i = 1, nnp%n_ele
1662 : rad = 0
1663 108 : ang = 0
1664 108 : DO j = 1, nnp%n_ele
1665 73 : funccut = -1.0_dp
1666 1246 : DO s = 1, nnp%n_rad(i)
1667 1211 : IF (nnp%rad(i)%ele(s) == nnp%ele(j)) THEN
1668 550 : IF (ABS(nnp%rad(i)%funccut(s) - funccut) > cutoff_eq_tol) THEN
1669 71 : rad = rad + 1
1670 71 : funccut = nnp%rad(i)%funccut(s)
1671 71 : nnp%rad(i)%symfgrp(rad)%cutoff = funccut
1672 71 : ALLOCATE (nnp%rad(i)%symfgrp(rad)%ele(1))
1673 71 : ALLOCATE (nnp%rad(i)%symfgrp(rad)%ele_ind(1))
1674 71 : nnp%rad(i)%symfgrp(rad)%ele(1) = nnp%ele(j)
1675 71 : nnp%rad(i)%symfgrp(rad)%ele_ind(1) = j
1676 : END IF
1677 550 : 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 125 : DO j = 1, nnp%n_ele
1682 222 : DO k = j, nnp%n_ele
1683 114 : funccut = -1.0_dp
1684 1507 : DO s = 1, nnp%n_ang(i)
1685 : IF ((nnp%ang(i)%ele1(s) == nnp%ele(j) .AND. &
1686 1320 : nnp%ang(i)%ele2(s) == nnp%ele(k)) .OR. &
1687 : (nnp%ang(i)%ele1(s) == nnp%ele(k) .AND. &
1688 114 : nnp%ang(i)%ele2(s) == nnp%ele(j))) THEN
1689 420 : IF (ABS(nnp%ang(i)%funccut(s) - funccut) > cutoff_eq_tol) THEN
1690 86 : ang = ang + 1
1691 86 : funccut = nnp%ang(i)%funccut(s)
1692 86 : nnp%ang(i)%symfgrp(ang)%cutoff = funccut
1693 86 : ALLOCATE (nnp%ang(i)%symfgrp(ang)%ele(2))
1694 86 : ALLOCATE (nnp%ang(i)%symfgrp(ang)%ele_ind(2))
1695 86 : nnp%ang(i)%symfgrp(ang)%ele(1) = nnp%ele(j)
1696 86 : nnp%ang(i)%symfgrp(ang)%ele(2) = nnp%ele(k)
1697 86 : nnp%ang(i)%symfgrp(ang)%ele_ind(1) = j
1698 86 : nnp%ang(i)%symfgrp(ang)%ele_ind(2) = k
1699 : END IF
1700 420 : 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 52 : DO i = 1, nnp%n_ele
1708 106 : DO j = 1, nnp%rad(i)%n_symfgrp
1709 213 : ALLOCATE (nnp%rad(i)%symfgrp(j)%symf(nnp%rad(i)%symfgrp(j)%n_symf))
1710 71 : rad = 0
1711 1223 : DO s = 1, nnp%n_rad(i)
1712 1188 : IF (nnp%rad(i)%ele(s) == nnp%rad(i)%symfgrp(j)%ele(1)) THEN
1713 550 : IF (ABS(nnp%rad(i)%funccut(s) - nnp%rad(i)%symfgrp(j)%cutoff) <= cutoff_eq_tol) THEN
1714 550 : rad = rad + 1
1715 550 : nnp%rad(i)%symfgrp(j)%symf(rad) = s
1716 : END IF
1717 : END IF
1718 : END DO
1719 : END DO
1720 138 : DO j = 1, nnp%ang(i)%n_symfgrp
1721 258 : ALLOCATE (nnp%ang(i)%symfgrp(j)%symf(nnp%ang(i)%symfgrp(j)%n_symf))
1722 86 : ang = 0
1723 1185 : DO s = 1, nnp%n_ang(i)
1724 : IF ((nnp%ang(i)%ele1(s) == nnp%ang(i)%symfgrp(j)%ele(1) .AND. &
1725 1064 : 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 86 : nnp%ang(i)%ele2(s) == nnp%ang(i)%symfgrp(j)%ele(1))) THEN
1728 420 : IF (ABS(nnp%ang(i)%funccut(s) - nnp%ang(i)%symfgrp(j)%cutoff) <= cutoff_eq_tol) THEN
1729 420 : ang = ang + 1
1730 420 : 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 52 : DO i = 1, nnp%n_ele
1745 138 : DO j = 1, nnp%ang(i)%n_symfgrp
1746 86 : n_symf = nnp%ang(i)%symfgrp(j)%n_symf
1747 258 : ALLOCATE (nnp%ang(i)%symfgrp(j)%pack_eta(n_symf))
1748 172 : ALLOCATE (nnp%ang(i)%symfgrp(j)%pack_zeta(n_symf))
1749 172 : ALLOCATE (nnp%ang(i)%symfgrp(j)%pack_lam(n_symf))
1750 172 : ALLOCATE (nnp%ang(i)%symfgrp(j)%pack_prefzeta(n_symf))
1751 258 : ALLOCATE (nnp%ang(i)%symfgrp(j)%pack_izeta(n_symf))
1752 172 : ALLOCATE (nnp%ang(i)%symfgrp(j)%pack_use_int_zeta(n_symf))
1753 541 : DO sf = 1, n_symf
1754 420 : m = nnp%ang(i)%symfgrp(j)%symf(sf)
1755 420 : eta_tmp = nnp%ang(i)%eta(m)
1756 420 : zeta_tmp = nnp%ang(i)%zeta(m)
1757 420 : nnp%ang(i)%symfgrp(j)%pack_eta(sf) = eta_tmp
1758 420 : nnp%ang(i)%symfgrp(j)%pack_zeta(sf) = zeta_tmp
1759 420 : nnp%ang(i)%symfgrp(j)%pack_lam(sf) = nnp%ang(i)%lam(m)
1760 420 : nnp%ang(i)%symfgrp(j)%pack_prefzeta(sf) = nnp%ang(i)%prefzeta(m)
1761 420 : izeta_tmp = NINT(zeta_tmp)
1762 420 : nnp%ang(i)%symfgrp(j)%pack_izeta(sf) = izeta_tmp
1763 : nnp%ang(i)%symfgrp(j)%pack_use_int_zeta(sf) = &
1764 506 : (REAL(izeta_tmp, dp) == zeta_tmp)
1765 : END DO
1766 : END DO
1767 : END DO
1768 :
1769 17 : 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 17 : 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 17 : NULLIFY (logger)
1791 17 : logger => cp_get_default_logger()
1792 :
1793 17 : my_label = TRIM(printtag)//"| "
1794 17 : IF (para_env%is_source()) THEN
1795 9 : unit_nr = cp_logger_get_default_unit_nr(logger)
1796 9 : WRITE (unit_nr, '(1X,A,1X,10(I2,1X))') TRIM(my_label)//" Activation functions:", nnp%actfnct(:)
1797 28 : DO i = 1, nnp%n_ele
1798 : WRITE (unit_nr, *) TRIM(my_label)//" short range atomic symmetry functions element "// &
1799 19 : nnp%ele(i)//":"
1800 313 : DO j = 1, nnp%n_rad(i)
1801 294 : 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 294 : nnp%rad(i)%ele(j), nnp%rad(i)%eta(j), &
1803 607 : nnp%rad(i)%rs(j), nnp%rad(i)%funccut(j)
1804 : END DO
1805 248 : 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 220 : TRIM(my_label), j, nnp%ele(i), 3, &
1808 220 : nnp%ang(i)%ele1(j), nnp%ang(i)%ele2(j), &
1809 220 : nnp%ang(i)%eta(j), nnp%ang(i)%lam(j), &
1810 459 : nnp%ang(i)%zeta(j), nnp%ang(i)%funccut(j)
1811 : END DO
1812 : END DO
1813 : END IF
1814 :
1815 17 : END SUBROUTINE nnp_write_acsf
1816 :
1817 : END MODULE nnp_acsf
|
