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