Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 : !> \date 2020-10-10
13 : ! **************************************************************************************************
14 : MODULE nnp_acsf
15 : USE cell_types, ONLY: cell_type,&
16 : pbc
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_environment_types, ONLY: nnp_acsf_ang_type,&
25 : nnp_acsf_rad_type,&
26 : nnp_cut_cos,&
27 : nnp_cut_tanh,&
28 : nnp_env_get,&
29 : nnp_neighbor_type,&
30 : nnp_type
31 : USE periodic_table, ONLY: get_ptable_info
32 : #include "./base/base_uses.f90"
33 :
34 : IMPLICIT NONE
35 :
36 : PRIVATE
37 :
38 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
39 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nnp_acsf'
40 : ! Public subroutines ***
41 : PUBLIC :: nnp_calc_acsf, &
42 : nnp_init_acsf_groups, &
43 : nnp_sort_acsf, &
44 : nnp_sort_ele, &
45 : nnp_write_acsf
46 :
47 : CONTAINS
48 :
49 : ! **************************************************************************************************
50 : !> \brief Calculate atom centered symmetry functions for given atom i
51 : !> \param nnp ...
52 : !> \param i ...
53 : !> \param dsymdxyz ...
54 : !> \param stress ...
55 : !> \date 2020-10-10
56 : !> \author Christoph Schran (christoph.schran@rub.de)
57 : ! **************************************************************************************************
58 248177 : SUBROUTINE nnp_calc_acsf(nnp, i, dsymdxyz, stress)
59 :
60 : TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
61 : INTEGER, INTENT(IN) :: i
62 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
63 : OPTIONAL :: dsymdxyz, stress
64 :
65 : CHARACTER(len=*), PARAMETER :: routineN = 'nnp_calc_acsf'
66 :
67 : INTEGER :: handle, handle_sf, ind, j, jj, k, kk, l, &
68 : m, off, s, sf
69 : REAL(KIND=dp) :: r1, r2, r3
70 248177 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: symtmp
71 248177 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: forcetmp
72 248177 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: force3tmp
73 : REAL(KIND=dp), DIMENSION(3) :: rvect1, rvect2, rvect3
74 992708 : TYPE(nnp_neighbor_type) :: neighbor
75 :
76 248177 : CALL timeset(routineN, handle)
77 :
78 : !determine index of atom type
79 248177 : ind = nnp%ele_ind(i)
80 :
81 : ! compute neighbors of atom i
82 248177 : CALL nnp_neighbor_create(nnp, ind, nnp%num_atoms, neighbor)
83 248177 : CALL nnp_compute_neighbors(nnp, neighbor, i)
84 :
85 : ! Reset y:
86 3722461 : nnp%rad(ind)%y = 0.0_dp
87 2067689 : nnp%ang(ind)%y = 0.0_dp
88 :
89 : !calc forces
90 248177 : IF (PRESENT(dsymdxyz)) THEN
91 : !loop over radial sym fnct grps
92 19105 : CALL timeset('nnp_acsf_radial', handle_sf)
93 60515 : DO s = 1, nnp%rad(ind)%n_symfgrp
94 124230 : ALLOCATE (symtmp(nnp%rad(ind)%symfgrp(s)%n_symf))
95 124230 : ALLOCATE (forcetmp(3, nnp%rad(ind)%symfgrp(s)%n_symf))
96 : !loop over associated neighbors
97 1401178 : DO j = 1, neighbor%n_rad(s)
98 5439072 : rvect1 = neighbor%dist_rad(1:3, j, s)
99 1359768 : r1 = neighbor%dist_rad(4, j, s)
100 1359768 : CALL nnp_calc_rad(nnp, ind, s, rvect1, r1, symtmp, forcetmp)
101 1359768 : jj = neighbor%ind_rad(j, s)
102 12256054 : DO sf = 1, nnp%rad(ind)%symfgrp(s)%n_symf
103 10854876 : m = nnp%rad(ind)%symfgrp(s)%symf(sf)
104 : ! update forces into dsymdxyz
105 43419504 : DO l = 1, 3
106 32564628 : dsymdxyz(l, m, i) = dsymdxyz(l, m, i) + forcetmp(l, sf)
107 43419504 : dsymdxyz(l, m, jj) = dsymdxyz(l, m, jj) - forcetmp(l, sf)
108 : END DO
109 10854876 : IF (PRESENT(stress)) THEN
110 7184768 : DO l = 1, 3
111 23350496 : stress(:, l, m) = stress(:, l, m) + rvect1(:)*forcetmp(l, sf)
112 : END DO
113 : END IF
114 12214644 : nnp%rad(ind)%y(m) = nnp%rad(ind)%y(m) + symtmp(sf)
115 : END DO
116 : END DO
117 41410 : DEALLOCATE (symtmp)
118 60515 : DEALLOCATE (forcetmp)
119 : END DO
120 19105 : CALL timestop(handle_sf)
121 : !loop over angular sym fnct grps
122 19105 : CALL timeset('nnp_acsf_angular', handle_sf)
123 19105 : off = nnp%n_rad(ind)
124 61550 : DO s = 1, nnp%ang(ind)%n_symfgrp
125 127335 : ALLOCATE (symtmp(nnp%ang(ind)%symfgrp(s)%n_symf))
126 127335 : ALLOCATE (force3tmp(3, 3, nnp%ang(ind)%symfgrp(s)%n_symf))
127 : !loop over associated neighbors
128 42445 : IF (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2)) THEN
129 774360 : DO j = 1, neighbor%n_ang1(s)
130 3016880 : rvect1 = neighbor%dist_ang1(1:3, j, s)
131 754220 : r1 = neighbor%dist_ang1(4, j, s)
132 19195970 : DO k = j + 1, neighbor%n_ang1(s)
133 73686440 : rvect2 = neighbor%dist_ang1(1:3, k, s)
134 18421610 : r2 = neighbor%dist_ang1(4, k, s)
135 73686440 : rvect3 = rvect2 - rvect1
136 73686440 : r3 = NORM2(rvect3(:))
137 19175830 : IF (r3 < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
138 8481295 : jj = neighbor%ind_ang1(j, s)
139 8481295 : kk = neighbor%ind_ang1(k, s)
140 : CALL nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, &
141 8481295 : r1, r2, r3, symtmp, force3tmp)
142 52399787 : DO sf = 1, nnp%ang(ind)%symfgrp(s)%n_symf
143 43918492 : m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
144 : ! update forces into dsymdxy
145 175673968 : DO l = 1, 3
146 : dsymdxyz(l, m, i) = dsymdxyz(l, m, i) &
147 131755476 : + force3tmp(l, 1, sf)
148 : dsymdxyz(l, m, jj) = dsymdxyz(l, m, jj) &
149 131755476 : + force3tmp(l, 2, sf)
150 : dsymdxyz(l, m, kk) = dsymdxyz(l, m, kk) &
151 175673968 : + force3tmp(l, 3, sf)
152 : END DO
153 43918492 : IF (PRESENT(stress)) THEN
154 29217640 : DO l = 1, 3
155 87652920 : stress(:, l, m) = stress(:, l, m) - rvect1(:)*force3tmp(l, 2, sf)
156 94957330 : stress(:, l, m) = stress(:, l, m) - rvect2(:)*force3tmp(l, 3, sf)
157 : END DO
158 : END IF
159 52399787 : nnp%ang(ind)%y(m - off) = nnp%ang(ind)%y(m - off) + symtmp(sf)
160 : END DO
161 : END IF
162 : END DO
163 : END DO
164 : ELSE
165 932948 : DO j = 1, neighbor%n_ang1(s)
166 3642572 : rvect1 = neighbor%dist_ang1(1:3, j, s)
167 910643 : r1 = neighbor%dist_ang1(4, j, s)
168 32831058 : DO k = 1, neighbor%n_ang2(s)
169 127592440 : rvect2 = neighbor%dist_ang2(1:3, k, s)
170 31898110 : r2 = neighbor%dist_ang2(4, k, s)
171 127592440 : rvect3 = rvect2 - rvect1
172 127592440 : r3 = NORM2(rvect3(:))
173 32808753 : IF (r3 < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
174 14833666 : jj = neighbor%ind_ang1(j, s)
175 14833666 : kk = neighbor%ind_ang1(k, s)
176 : CALL nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, &
177 14833666 : r1, r2, r3, symtmp, force3tmp)
178 : !loop over associated sym fncts
179 104150027 : DO sf = 1, nnp%ang(ind)%symfgrp(s)%n_symf
180 89316361 : m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
181 89316361 : jj = neighbor%ind_ang1(j, s)
182 89316361 : kk = neighbor%ind_ang2(k, s)
183 : ! update forces into dsymdxy
184 357265444 : DO l = 1, 3
185 : dsymdxyz(l, m, i) = dsymdxyz(l, m, i) &
186 267949083 : + force3tmp(l, 1, sf)
187 : dsymdxyz(l, m, jj) = dsymdxyz(l, m, jj) &
188 267949083 : + force3tmp(l, 2, sf)
189 : dsymdxyz(l, m, kk) = dsymdxyz(l, m, kk) &
190 357265444 : + force3tmp(l, 3, sf)
191 : END DO
192 89316361 : IF (PRESENT(stress)) THEN
193 59415432 : DO l = 1, 3
194 178246296 : stress(:, l, m) = stress(:, l, m) - rvect1(:)*force3tmp(l, 2, sf)
195 193100154 : stress(:, l, m) = stress(:, l, m) - rvect2(:)*force3tmp(l, 3, sf)
196 : END DO
197 : END IF
198 104150027 : nnp%ang(ind)%y(m - off) = nnp%ang(ind)%y(m - off) + symtmp(sf)
199 : END DO
200 : END IF
201 : END DO
202 : END DO
203 : END IF
204 42445 : DEALLOCATE (symtmp)
205 61550 : DEALLOCATE (force3tmp)
206 : END DO
207 19105 : CALL timestop(handle_sf)
208 : !no forces
209 : ELSE
210 : !loop over radial sym fnct grps
211 229072 : CALL timeset('nnp_acsf_radial', handle_sf)
212 794360 : DO s = 1, nnp%rad(ind)%n_symfgrp
213 1695864 : ALLOCATE (symtmp(nnp%rad(ind)%symfgrp(s)%n_symf))
214 : !loop over associated neighbors
215 2495804 : DO j = 1, neighbor%n_rad(s)
216 7722064 : rvect1 = neighbor%dist_rad(1:3, j, s)
217 1930516 : r1 = neighbor%dist_rad(4, j, s)
218 1930516 : CALL nnp_calc_rad(nnp, ind, s, rvect1, r1, symtmp)
219 17187322 : DO sf = 1, nnp%rad(ind)%symfgrp(s)%n_symf
220 14691518 : m = nnp%rad(ind)%symfgrp(s)%symf(sf)
221 16622034 : nnp%rad(ind)%y(m) = nnp%rad(ind)%y(m) + symtmp(sf)
222 : END DO
223 : END DO
224 794360 : DEALLOCATE (symtmp)
225 : END DO
226 229072 : CALL timestop(handle_sf)
227 : !loop over angular sym fnct grps
228 229072 : CALL timeset('nnp_acsf_angular', handle_sf)
229 229072 : off = nnp%n_rad(ind)
230 692144 : DO s = 1, nnp%ang(ind)%n_symfgrp
231 1389216 : ALLOCATE (symtmp(nnp%ang(ind)%symfgrp(s)%n_symf))
232 : !loop over associated neighbors
233 463072 : IF (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2)) THEN
234 1121116 : DO j = 1, neighbor%n_ang1(s)
235 3977040 : rvect1 = neighbor%dist_ang1(1:3, j, s)
236 994260 : r1 = neighbor%dist_ang1(4, j, s)
237 22659329 : DO k = j + 1, neighbor%n_ang1(s)
238 86152852 : rvect2 = neighbor%dist_ang1(1:3, k, s)
239 21538213 : r2 = neighbor%dist_ang1(4, k, s)
240 86152852 : rvect3 = rvect2 - rvect1
241 86152852 : r3 = NORM2(rvect3(:))
242 22532473 : IF (r3 < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
243 9944571 : CALL nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, r1, r2, r3, symtmp)
244 61315369 : DO sf = 1, nnp%ang(ind)%symfgrp(s)%n_symf
245 51370798 : m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
246 61315369 : nnp%ang(ind)%y(m - off) = nnp%ang(ind)%y(m - off) + symtmp(sf)
247 : END DO
248 : END IF
249 : END DO
250 : END DO
251 : ELSE
252 1719216 : DO j = 1, neighbor%n_ang1(s)
253 5532000 : rvect1 = neighbor%dist_ang1(1:3, j, s)
254 1383000 : r1 = neighbor%dist_ang1(4, j, s)
255 39030025 : DO k = 1, neighbor%n_ang2(s)
256 149243236 : rvect2 = neighbor%dist_ang2(1:3, k, s)
257 37310809 : r2 = neighbor%dist_ang2(4, k, s)
258 149243236 : rvect3 = rvect2 - rvect1
259 149243236 : r3 = NORM2(rvect3(:))
260 38693809 : IF (r3 < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
261 17429498 : CALL nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, r1, r2, r3, symtmp)
262 : !loop over associated sym fncts
263 121983157 : DO sf = 1, nnp%ang(ind)%symfgrp(s)%n_symf
264 104553659 : m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
265 121983157 : nnp%ang(ind)%y(m - off) = nnp%ang(ind)%y(m - off) + symtmp(sf)
266 : END DO
267 : END IF
268 : END DO
269 : END DO
270 : END IF
271 692144 : DEALLOCATE (symtmp)
272 : END DO
273 229072 : CALL timestop(handle_sf)
274 : END IF
275 :
276 : !check extrapolation
277 248177 : CALL nnp_check_extrapolation(nnp, ind)
278 :
279 248177 : IF (PRESENT(dsymdxyz)) THEN
280 19105 : IF (PRESENT(stress)) THEN
281 2112 : CALL nnp_scale_acsf(nnp, ind, dsymdxyz, stress)
282 : ELSE
283 16993 : CALL nnp_scale_acsf(nnp, ind, dsymdxyz)
284 : END IF
285 : ELSE
286 229072 : CALL nnp_scale_acsf(nnp, ind)
287 : END IF
288 :
289 248177 : CALL nnp_neighbor_release(neighbor)
290 248177 : CALL timestop(handle)
291 :
292 248177 : END SUBROUTINE nnp_calc_acsf
293 :
294 : ! **************************************************************************************************
295 : !> \brief Check if the nnp is extrapolating
296 : !> \param nnp ...
297 : !> \param ind ...
298 : !> \date 2020-10-10
299 : !> \author Christoph Schran (christoph.schran@rub.de)
300 : ! **************************************************************************************************
301 248177 : SUBROUTINE nnp_check_extrapolation(nnp, ind)
302 :
303 : TYPE(nnp_type), INTENT(INOUT) :: nnp
304 : INTEGER, INTENT(IN) :: ind
305 :
306 : REAL(KIND=dp), PARAMETER :: threshold = 0.0001_dp
307 :
308 : INTEGER :: j
309 : LOGICAL :: extrapolate
310 :
311 248177 : extrapolate = nnp%output_expol
312 :
313 3722461 : DO j = 1, nnp%n_rad(ind)
314 3474284 : IF (nnp%rad(ind)%y(j) - &
315 248177 : nnp%rad(ind)%loc_max(j) > threshold) THEN
316 : extrapolate = .TRUE.
317 3474284 : ELSE IF (-nnp%rad(ind)%y(j) + &
318 : nnp%rad(ind)%loc_min(j) > threshold) THEN
319 164 : extrapolate = .TRUE.
320 : END IF
321 : END DO
322 2067689 : DO j = 1, nnp%n_ang(ind)
323 1819512 : IF (nnp%ang(ind)%y(j) - &
324 248177 : nnp%ang(ind)%loc_max(j) > threshold) THEN
325 : extrapolate = .TRUE.
326 1819512 : ELSE IF (-nnp%ang(ind)%y(j) + &
327 : nnp%ang(ind)%loc_min(j) > threshold) THEN
328 30 : extrapolate = .TRUE.
329 : END IF
330 : END DO
331 :
332 248177 : nnp%output_expol = extrapolate
333 :
334 248177 : END SUBROUTINE nnp_check_extrapolation
335 :
336 : ! **************************************************************************************************
337 : !> \brief Scale and center symetry functions (and gradients)
338 : !> \param nnp ...
339 : !> \param ind ...
340 : !> \param dsymdxyz ...
341 : !> \param stress ...
342 : !> \date 2020-10-10
343 : !> \author Christoph Schran (christoph.schran@rub.de)
344 : ! **************************************************************************************************
345 248177 : SUBROUTINE nnp_scale_acsf(nnp, ind, dsymdxyz, stress)
346 :
347 : TYPE(nnp_type), INTENT(INOUT) :: nnp
348 : INTEGER, INTENT(IN) :: ind
349 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT), &
350 : OPTIONAL :: dsymdxyz, stress
351 :
352 : INTEGER :: j, k, off
353 :
354 248177 : IF (nnp%center_acsf) THEN
355 3722461 : DO j = 1, nnp%n_rad(ind)
356 : nnp%arc(ind)%layer(1)%node(j) = &
357 3722461 : (nnp%rad(ind)%y(j) - nnp%rad(ind)%loc_av(j))
358 : END DO
359 248177 : off = nnp%n_rad(ind)
360 2067689 : DO j = 1, nnp%n_ang(ind)
361 : nnp%arc(ind)%layer(1)%node(j + off) = &
362 2067689 : (nnp%ang(ind)%y(j) - nnp%ang(ind)%loc_av(j))
363 : END DO
364 :
365 248177 : IF (nnp%scale_acsf) THEN
366 3722461 : DO j = 1, nnp%n_rad(ind)
367 : nnp%arc(ind)%layer(1)%node(j) = &
368 : nnp%arc(ind)%layer(1)%node(j)/ &
369 : (nnp%rad(ind)%loc_max(j) - nnp%rad(ind)%loc_min(j))* &
370 3722461 : (nnp%scmax - nnp%scmin) + nnp%scmin
371 : END DO
372 248177 : off = nnp%n_rad(ind)
373 2067689 : DO j = 1, nnp%n_ang(ind)
374 : nnp%arc(ind)%layer(1)%node(j + off) = &
375 : nnp%arc(ind)%layer(1)%node(j + off)/ &
376 : (nnp%ang(ind)%loc_max(j) - nnp%ang(ind)%loc_min(j))* &
377 2067689 : (nnp%scmax - nnp%scmin) + nnp%scmin
378 : END DO
379 : END IF
380 0 : ELSE IF (nnp%scale_acsf) THEN
381 0 : DO j = 1, nnp%n_rad(ind)
382 : nnp%arc(ind)%layer(1)%node(j) = &
383 : (nnp%rad(ind)%y(j) - nnp%rad(ind)%loc_min(j))/ &
384 : (nnp%rad(ind)%loc_max(j) - nnp%rad(ind)%loc_min(j))* &
385 0 : (nnp%scmax - nnp%scmin) + nnp%scmin
386 : END DO
387 0 : off = nnp%n_rad(ind)
388 0 : DO j = 1, nnp%n_ang(ind)
389 : nnp%arc(ind)%layer(1)%node(j + off) = &
390 : (nnp%ang(ind)%y(j) - nnp%ang(ind)%loc_min(j))/ &
391 : (nnp%ang(ind)%loc_max(j) - nnp%ang(ind)%loc_min(j))* &
392 0 : (nnp%scmax - nnp%scmin) + nnp%scmin
393 : END DO
394 0 : ELSE IF (nnp%scale_sigma_acsf) THEN
395 0 : DO j = 1, nnp%n_rad(ind)
396 : nnp%arc(ind)%layer(1)%node(j) = &
397 : (nnp%rad(ind)%y(j) - nnp%rad(ind)%loc_av(j))/ &
398 : nnp%rad(ind)%sigma(j)* &
399 0 : (nnp%scmax - nnp%scmin) + nnp%scmin
400 : END DO
401 0 : off = nnp%n_rad(ind)
402 0 : DO j = 1, nnp%n_ang(ind)
403 : nnp%arc(ind)%layer(1)%node(j + off) = &
404 : (nnp%ang(ind)%y(j) - nnp%ang(ind)%loc_av(j))/ &
405 : nnp%ang(ind)%sigma(j)* &
406 0 : (nnp%scmax - nnp%scmin) + nnp%scmin
407 : END DO
408 : ELSE
409 0 : DO j = 1, nnp%n_rad(ind)
410 0 : nnp%arc(ind)%layer(1)%node(j) = nnp%rad(ind)%y(j)
411 : END DO
412 0 : off = nnp%n_rad(ind)
413 0 : DO j = 1, nnp%n_ang(ind)
414 0 : nnp%arc(ind)%layer(1)%node(j + off) = nnp%ang(ind)%y(j)
415 : END DO
416 : END IF
417 :
418 248177 : IF (PRESENT(dsymdxyz)) THEN
419 19105 : IF (nnp%scale_acsf) THEN
420 2477828 : DO k = 1, nnp%num_atoms
421 41759796 : DO j = 1, nnp%n_rad(ind)
422 : dsymdxyz(:, j, k) = dsymdxyz(:, j, k)/ &
423 : (nnp%rad(ind)%loc_max(j) - &
424 : nnp%rad(ind)%loc_min(j))* &
425 159586595 : (nnp%scmax - nnp%scmin)
426 : END DO
427 : END DO
428 19105 : off = nnp%n_rad(ind)
429 2477828 : DO k = 1, nnp%num_atoms
430 31848104 : DO j = 1, nnp%n_ang(ind)
431 : dsymdxyz(:, j + off, k) = dsymdxyz(:, j + off, k)/ &
432 : (nnp%ang(ind)%loc_max(j) - &
433 : nnp%ang(ind)%loc_min(j))* &
434 119939827 : (nnp%scmax - nnp%scmin)
435 : END DO
436 : END DO
437 0 : ELSE IF (nnp%scale_sigma_acsf) THEN
438 0 : DO k = 1, nnp%num_atoms
439 0 : DO j = 1, nnp%n_rad(ind)
440 : dsymdxyz(:, j, k) = dsymdxyz(:, j, k)/ &
441 : nnp%rad(ind)%sigma(j)* &
442 0 : (nnp%scmax - nnp%scmin)
443 : END DO
444 : END DO
445 0 : off = nnp%n_rad(ind)
446 0 : DO k = 1, nnp%num_atoms
447 0 : DO j = 1, nnp%n_ang(ind)
448 : dsymdxyz(:, j + off, k) = dsymdxyz(:, j + off, k)/ &
449 : nnp%ang(ind)%sigma(j)* &
450 0 : (nnp%scmax - nnp%scmin)
451 : END DO
452 : END DO
453 : END IF
454 : END IF
455 :
456 248177 : IF (PRESENT(stress)) THEN
457 2112 : IF (nnp%scale_acsf) THEN
458 35904 : DO j = 1, nnp%n_rad(ind)
459 : stress(:, :, j) = stress(:, :, j)/ &
460 : (nnp%rad(ind)%loc_max(j) - &
461 : nnp%rad(ind)%loc_min(j))* &
462 441408 : (nnp%scmax - nnp%scmin)
463 : END DO
464 2112 : off = nnp%n_rad(ind)
465 27456 : DO j = 1, nnp%n_ang(ind)
466 : stress(:, :, j + off) = stress(:, :, j + off)/ &
467 : (nnp%ang(ind)%loc_max(j) - &
468 : nnp%ang(ind)%loc_min(j))* &
469 331584 : (nnp%scmax - nnp%scmin)
470 : END DO
471 0 : ELSE IF (nnp%scale_sigma_acsf) THEN
472 0 : DO j = 1, nnp%n_rad(ind)
473 : stress(:, :, j) = stress(:, :, j)/ &
474 : nnp%rad(ind)%sigma(j)* &
475 0 : (nnp%scmax - nnp%scmin)
476 : END DO
477 0 : off = nnp%n_rad(ind)
478 0 : DO j = 1, nnp%n_ang(ind)
479 : stress(:, :, j + off) = stress(:, :, j + off)/ &
480 : nnp%ang(ind)%sigma(j)* &
481 0 : (nnp%scmax - nnp%scmin)
482 : END DO
483 : END IF
484 : END IF
485 :
486 248177 : END SUBROUTINE nnp_scale_acsf
487 :
488 : ! **************************************************************************************************
489 : !> \brief Calculate radial symmetry function and gradient (optinal)
490 : !> for given displacment vecotr rvect of atom i and j
491 : !> \param nnp ...
492 : !> \param ind ...
493 : !> \param s ...
494 : !> \param rvect ...
495 : !> \param r ...
496 : !> \param sym ...
497 : !> \param force ...
498 : !> \date 2020-10-10
499 : !> \author Christoph Schran (christoph.schran@rub.de)
500 : ! **************************************************************************************************
501 3290284 : SUBROUTINE nnp_calc_rad(nnp, ind, s, rvect, r, sym, force)
502 :
503 : TYPE(nnp_type), INTENT(IN) :: nnp
504 : INTEGER, INTENT(IN) :: ind, s
505 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rvect
506 : REAL(KIND=dp), INTENT(IN) :: r
507 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: sym
508 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), &
509 : OPTIONAL :: force
510 :
511 : INTEGER :: k, sf
512 : REAL(KIND=dp) :: dfcutdr, dsymdr, eta, fcut, funccut, rs, &
513 : tmp
514 : REAL(KIND=dp), DIMENSION(3) :: drdx
515 :
516 : !init
517 3290284 : drdx = 0.0_dp
518 3290284 : fcut = 0.0_dp
519 3290284 : dfcutdr = 0.0_dp
520 :
521 : !Calculate cutoff function and partial derivative
522 3290284 : funccut = nnp%rad(ind)%symfgrp(s)%cutoff !cutoff
523 3659254 : SELECT CASE (nnp%cut_type)
524 : CASE (nnp_cut_cos)
525 368970 : tmp = pi*r/funccut
526 368970 : fcut = 0.5_dp*(COS(tmp) + 1.0_dp)
527 368970 : IF (PRESENT(force)) THEN
528 10956 : dfcutdr = 0.5_dp*(-SIN(tmp))*(pi/funccut)
529 : END IF
530 : CASE (nnp_cut_tanh)
531 2921314 : tmp = TANH(1.0_dp - r/funccut)
532 2921314 : fcut = tmp**3
533 2921314 : IF (PRESENT(force)) THEN
534 1348812 : dfcutdr = (-3.0_dp/funccut)*(tmp**2 - tmp**4)
535 : END IF
536 : CASE DEFAULT
537 3290284 : CPABORT("NNP| Cutoff function unknown")
538 : END SELECT
539 :
540 7369588 : IF (PRESENT(force)) drdx(:) = rvect(:)/r
541 :
542 : !loop over sym fncts of sym fnct group s
543 28836678 : DO sf = 1, nnp%rad(ind)%symfgrp(s)%n_symf
544 25546394 : k = nnp%rad(ind)%symfgrp(s)%symf(sf) !symf indice
545 25546394 : eta = nnp%rad(ind)%eta(k) !eta
546 25546394 : rs = nnp%rad(ind)%rs(k) !rshift
547 :
548 : ! Calculate radial symmetry function
549 25546394 : sym(sf) = EXP(-eta*(r - rs)**2)
550 :
551 : ! Calculate partial derivatives of symmetry function and distance
552 : ! and combine them to obtain force
553 25546394 : IF (PRESENT(force)) THEN
554 10854876 : dsymdr = sym(sf)*(-2.0_dp*eta*(r - rs))
555 43419504 : force(:, sf) = fcut*dsymdr*drdx(:) + sym(sf)*dfcutdr*drdx(:)
556 : END IF
557 :
558 : ! combine radial symmetry function and cutoff function
559 28836678 : sym(sf) = sym(sf)*fcut
560 : END DO
561 :
562 3290284 : END SUBROUTINE nnp_calc_rad
563 :
564 : ! **************************************************************************************************
565 : !> \brief Calculate angular symmetry function and gradient (optinal)
566 : !> for given displacment vectors rvect1,2,3 of atom i,j and k
567 : !> \param nnp ...
568 : !> \param ind ...
569 : !> \param s ...
570 : !> \param rvect1 ...
571 : !> \param rvect2 ...
572 : !> \param rvect3 ...
573 : !> \param r1 ...
574 : !> \param r2 ...
575 : !> \param r3 ...
576 : !> \param sym ...
577 : !> \param force ...
578 : !> \date 2020-10-10
579 : !> \author Christoph Schran (christoph.schran@rub.de)
580 : ! **************************************************************************************************
581 50689030 : SUBROUTINE nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, r1, r2, r3, sym, force)
582 :
583 : TYPE(nnp_type), INTENT(IN) :: nnp
584 : INTEGER, INTENT(IN) :: ind, s
585 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rvect1, rvect2, rvect3
586 : REAL(KIND=dp), INTENT(IN) :: r1, r2, r3
587 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: sym
588 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT), &
589 : OPTIONAL :: force
590 :
591 : INTEGER :: i, m, sf
592 : REAL(KIND=dp) :: angular, costheta, dfcutdr1, dfcutdr2, dfcutdr3, dsymdr1, dsymdr2, dsymdr3, &
593 : eta, f, fcut, fcut1, fcut2, fcut3, ftot, g, lam, pref, prefzeta, rsqr1, rsqr2, rsqr3, &
594 : symtmp, tmp, tmp1, tmp2, tmp3, tmpzeta, zeta
595 : REAL(KIND=dp), DIMENSION(3) :: dangulardx1, dangulardx2, dangulardx3, dcosthetadx1, &
596 : dcosthetadx2, dcosthetadx3, dfdx1, dfdx2, dfdx3, dgdx1, dgdx2, dgdx3, dr1dx, dr2dx, dr3dx
597 :
598 50689030 : rsqr1 = r1**2
599 50689030 : rsqr2 = r2**2
600 50689030 : rsqr3 = r3**2
601 :
602 : !init
603 50689030 : dangulardx1 = 0.0_dp
604 50689030 : dangulardx2 = 0.0_dp
605 50689030 : dangulardx3 = 0.0_dp
606 50689030 : dr1dx = 0.0_dp
607 50689030 : dr2dx = 0.0_dp
608 50689030 : dr3dx = 0.0_dp
609 50689030 : ftot = 0.0_dp
610 50689030 : dfcutdr1 = 0.0_dp
611 50689030 : dfcutdr2 = 0.0_dp
612 50689030 : dfcutdr3 = 0.0_dp
613 :
614 : ! Calculate cos(theta)
615 : ! use law of cosine for theta
616 : ! cos(ang(r1,r2)) = (r3**2 - r1**2 - r2**2)/(-2*r1*r2)
617 : ! | f | g |
618 50689030 : f = (rsqr3 - rsqr1 - rsqr2)
619 50689030 : g = -2.0_dp*r1*r2
620 50689030 : costheta = f/g
621 :
622 : ! Calculate cutoff function and partial derivatives
623 50689030 : fcut = nnp%ang(ind)%symfgrp(s)%cutoff !cutoff
624 50894834 : SELECT CASE (nnp%cut_type)
625 : CASE (nnp_cut_cos)
626 205804 : tmp1 = pi*r1/fcut
627 205804 : tmp2 = pi*r2/fcut
628 205804 : tmp3 = pi*r3/fcut
629 205804 : fcut1 = 0.5_dp*(COS(tmp1) + 1.0_dp)
630 205804 : fcut2 = 0.5_dp*(COS(tmp2) + 1.0_dp)
631 205804 : fcut3 = 0.5_dp*(COS(tmp3) + 1.0_dp)
632 205804 : ftot = fcut1*fcut2*fcut3
633 205804 : IF (PRESENT(force)) THEN
634 6242 : pref = 0.5_dp*(pi/fcut)
635 6242 : dfcutdr1 = pref*(-SIN(tmp1))*fcut2*fcut3
636 6242 : dfcutdr2 = pref*(-SIN(tmp2))*fcut1*fcut3
637 6242 : dfcutdr3 = pref*(-SIN(tmp3))*fcut1*fcut2
638 : END IF
639 : CASE (nnp_cut_tanh)
640 50483226 : tmp1 = TANH(1.0_dp - r1/fcut)
641 50483226 : tmp2 = TANH(1.0_dp - r2/fcut)
642 50483226 : tmp3 = TANH(1.0_dp - r3/fcut)
643 50483226 : fcut1 = tmp1**3
644 50483226 : fcut2 = tmp2**3
645 50483226 : fcut3 = tmp3**3
646 50483226 : ftot = fcut1*fcut2*fcut3
647 50483226 : IF (PRESENT(force)) THEN
648 23308719 : pref = -3.0_dp/fcut
649 23308719 : dfcutdr1 = pref*(tmp1**2 - tmp1**4)*fcut2*fcut3
650 23308719 : dfcutdr2 = pref*(tmp2**2 - tmp2**4)*fcut1*fcut3
651 23308719 : dfcutdr3 = pref*(tmp3**2 - tmp3**4)*fcut1*fcut2
652 : END IF
653 : CASE DEFAULT
654 50689030 : CPABORT("NNP| Cutoff function unknown")
655 : END SELECT
656 :
657 50689030 : IF (PRESENT(force)) THEN
658 93259844 : dr1dx(:) = rvect1(:)/r1
659 93259844 : dr2dx(:) = rvect2(:)/r2
660 93259844 : dr3dx(:) = rvect3(:)/r3
661 : END IF
662 :
663 : !loop over associated sym fncts
664 339848340 : DO sf = 1, nnp%ang(ind)%symfgrp(s)%n_symf
665 289159310 : m = nnp%ang(ind)%symfgrp(s)%symf(sf)
666 289159310 : lam = nnp%ang(ind)%lam(m) !lambda
667 289159310 : zeta = nnp%ang(ind)%zeta(m) !zeta
668 289159310 : prefzeta = nnp%ang(ind)%prefzeta(m) ! 2**(1-zeta)
669 289159310 : eta = nnp%ang(ind)%eta(m) !eta
670 :
671 289159310 : tmp = (1.0_dp + lam*costheta)
672 :
673 289159310 : IF (tmp <= 0.0_dp) THEN
674 0 : sym(sf) = 0.0_dp
675 0 : IF (PRESENT(force)) force(:, :, sf) = 0.0_dp
676 : CYCLE
677 : END IF
678 :
679 : ! Calculate symmetry function
680 : ! Calculate angular symmetry function
681 : ! ang = (1+lam*cos(theta_ijk))**zeta
682 289159310 : i = NINT(zeta)
683 289159310 : IF (1.0_dp*i == zeta) THEN
684 289159310 : tmpzeta = tmp**(i - 1) ! integer power is a LOT faster
685 : ELSE
686 0 : tmpzeta = tmp**(zeta - 1.0_dp)
687 : END IF
688 289159310 : angular = tmpzeta*tmp
689 : ! exponential part:
690 : ! exp(-eta*(r1**2+r2**2+r3**2))
691 289159310 : symtmp = EXP(-eta*(rsqr1 + rsqr2 + rsqr3))
692 :
693 : ! Partial derivatives (f/g)' = (f'g - fg')/g^2
694 289159310 : IF (PRESENT(force)) THEN
695 133234853 : pref = zeta*(tmpzeta)/g**2
696 532939412 : DO i = 1, 3
697 399704559 : dfdx1(i) = -2.0_dp*lam*(rvect1(i) + rvect2(i))
698 399704559 : dfdx2(i) = 2.0_dp*lam*(rvect3(i) + rvect1(i))
699 399704559 : dfdx3(i) = 2.0_dp*lam*(rvect2(i) - rvect3(i))
700 :
701 399704559 : tmp1 = 2.0_dp*r2*dr1dx(i)
702 399704559 : tmp2 = 2.0_dp*r1*dr2dx(i)
703 399704559 : dgdx1(i) = -(tmp1 + tmp2)
704 : dgdx2(i) = tmp1
705 : dgdx3(i) = tmp2
706 :
707 399704559 : dcosthetadx1(i) = dfdx1(i)*g - lam*f*dgdx1(i)
708 399704559 : dcosthetadx2(i) = dfdx2(i)*g - lam*f*dgdx2(i)
709 399704559 : dcosthetadx3(i) = dfdx3(i)*g - lam*f*dgdx3(i)
710 :
711 399704559 : dangulardx1(i) = pref*dcosthetadx1(i)
712 399704559 : dangulardx2(i) = pref*dcosthetadx2(i)
713 532939412 : dangulardx3(i) = pref*dcosthetadx3(i)
714 : END DO
715 :
716 : ! Calculate partial derivatives of exponential part and distance
717 : ! and combine partial derivatives to obtain force
718 133234853 : pref = prefzeta
719 133234853 : tmp = -2.0_dp*symtmp*eta
720 133234853 : dsymdr1 = tmp*r1
721 133234853 : dsymdr2 = tmp*r2
722 133234853 : dsymdr3 = tmp*r3
723 :
724 : ! G(r1,r2,r3) = pref*angular(r1,r2,r3)*exp(r1,r2,r3)*fcut(r1,r2,r3)
725 : ! dG/dx1 = (dangular/dx1* exp * fcut +
726 : ! angular * dexp/dx1* fcut +
727 : ! angular * exp * dfcut/dx1)*pref
728 : ! dr1/dx1 = -dr1/dx2
729 133234853 : tmp = pref*symtmp*ftot
730 133234853 : tmp1 = pref*angular*(ftot*dsymdr1 + dfcutdr1*symtmp)
731 133234853 : tmp2 = pref*angular*(ftot*dsymdr2 + dfcutdr2*symtmp)
732 133234853 : tmp3 = pref*angular*(ftot*dsymdr3 + dfcutdr3*symtmp)
733 532939412 : DO i = 1, 3
734 399704559 : force(i, 1, sf) = tmp*dangulardx1(i) + tmp1*dr1dx(i) + tmp2*dr2dx(i)
735 399704559 : force(i, 2, sf) = tmp*dangulardx2(i) - tmp1*dr1dx(i) + tmp3*dr3dx(i)
736 532939412 : force(i, 3, sf) = tmp*dangulardx3(i) - tmp2*dr2dx(i) - tmp3*dr3dx(i)
737 : END DO
738 : END IF
739 :
740 : ! Don't forget prefactor: 2**(1-ang%zeta)
741 289159310 : pref = prefzeta
742 : ! combine angular and exponential part with cutoff function
743 339848340 : sym(sf) = pref*angular*symtmp*ftot
744 : END DO
745 :
746 50689030 : END SUBROUTINE nnp_calc_ang
747 :
748 : ! **************************************************************************************************
749 : !> \brief Sort element array according to atomic number
750 : !> \param ele ...
751 : !> \param nuc_ele ...
752 : !> \date 2020-10-10
753 : !> \author Christoph Schran (christoph.schran@rub.de)
754 : ! **************************************************************************************************
755 15 : SUBROUTINE nnp_sort_ele(ele, nuc_ele)
756 : CHARACTER(len=2), DIMENSION(:), INTENT(INOUT) :: ele
757 : INTEGER, DIMENSION(:), INTENT(INOUT) :: nuc_ele
758 :
759 : CHARACTER(len=2) :: tmp_ele
760 : INTEGER :: i, j, loc, minimum, tmp_nuc_ele
761 :
762 : ! Determine atomic number
763 46 : DO i = 1, SIZE(ele)
764 46 : CALL get_ptable_info(ele(i), number=nuc_ele(i))
765 : END DO
766 :
767 : ! Sort both arrays
768 31 : DO i = 1, SIZE(ele) - 1
769 16 : minimum = nuc_ele(i)
770 16 : loc = i
771 33 : DO j = i + 1, SIZE(ele)
772 33 : IF (nuc_ele(j) < minimum) THEN
773 16 : loc = j
774 16 : minimum = nuc_ele(j)
775 : END IF
776 : END DO
777 16 : tmp_nuc_ele = nuc_ele(i)
778 16 : nuc_ele(i) = nuc_ele(loc)
779 16 : nuc_ele(loc) = tmp_nuc_ele
780 :
781 16 : tmp_ele = ele(i)
782 16 : ele(i) = ele(loc)
783 31 : ele(loc) = tmp_ele
784 :
785 : END DO
786 :
787 15 : END SUBROUTINE nnp_sort_ele
788 :
789 : ! **************************************************************************************************
790 : !> \brief Sort symmetry functions according to different criteria
791 : !> \param nnp ...
792 : !> \date 2020-10-10
793 : !> \author Christoph Schran (christoph.schran@rub.de)
794 : ! **************************************************************************************************
795 15 : SUBROUTINE nnp_sort_acsf(nnp)
796 : TYPE(nnp_type), INTENT(INOUT) :: nnp
797 :
798 : INTEGER :: i, j, k, loc
799 :
800 : ! First sort is according to symmetry function type
801 : ! This is done manually, since data structures are separate
802 : ! Note: Bubble sort is OK here, since small sort + special
803 46 : DO i = 1, nnp%n_ele
804 : ! sort by cutoff
805 : ! rad
806 486 : DO j = 1, nnp%n_rad(i) - 1
807 455 : loc = j
808 4051 : DO k = j + 1, nnp%n_rad(i)
809 4051 : IF (nnp%rad(i)%funccut(loc) > nnp%rad(i)%funccut(k)) THEN
810 6 : loc = k
811 : END IF
812 : END DO
813 : ! swap symfnct
814 486 : CALL nnp_swaprad(nnp%rad(i), j, loc)
815 : END DO
816 :
817 : ! sort by eta
818 : ! rad
819 486 : DO j = 1, nnp%n_rad(i) - 1
820 455 : loc = j
821 4051 : DO k = j + 1, nnp%n_rad(i)
822 3596 : IF (nnp%rad(i)%funccut(loc) == nnp%rad(i)%funccut(k) .AND. &
823 455 : nnp%rad(i)%eta(loc) > nnp%rad(i)%eta(k)) THEN
824 484 : loc = k
825 : END IF
826 : END DO
827 : ! swap symfnct
828 486 : CALL nnp_swaprad(nnp%rad(i), j, loc)
829 : END DO
830 :
831 : ! sort by rshift
832 : ! rad
833 486 : DO j = 1, nnp%n_rad(i) - 1
834 455 : loc = j
835 4051 : DO k = j + 1, nnp%n_rad(i)
836 : IF (nnp%rad(i)%funccut(loc) == nnp%rad(i)%funccut(k) .AND. &
837 3596 : nnp%rad(i)%eta(loc) == nnp%rad(i)%eta(k) .AND. &
838 455 : nnp%rad(i)%rs(loc) > nnp%rad(i)%rs(k)) THEN
839 56 : loc = k
840 : END IF
841 : END DO
842 : ! swap symfnct
843 486 : CALL nnp_swaprad(nnp%rad(i), j, loc)
844 : END DO
845 :
846 : ! sort by ele
847 : ! rad
848 486 : DO j = 1, nnp%n_rad(i) - 1
849 455 : loc = j
850 4051 : DO k = j + 1, nnp%n_rad(i)
851 : IF (nnp%rad(i)%funccut(loc) == nnp%rad(i)%funccut(k) .AND. &
852 : nnp%rad(i)%eta(loc) == nnp%rad(i)%eta(k) .AND. &
853 3596 : nnp%rad(i)%rs(loc) == nnp%rad(i)%rs(k) .AND. &
854 455 : nnp%rad(i)%nuc_ele(loc) > nnp%rad(i)%nuc_ele(k)) THEN
855 0 : loc = k
856 : END IF
857 : END DO
858 : ! swap symfnct
859 486 : CALL nnp_swaprad(nnp%rad(i), j, loc)
860 : END DO
861 :
862 : ! ang
863 : ! sort by cutoff
864 370 : DO j = 1, nnp%n_ang(i) - 1
865 339 : loc = j
866 2440 : DO k = j + 1, nnp%n_ang(i)
867 2440 : IF (nnp%ang(i)%funccut(loc) > nnp%ang(i)%funccut(k)) THEN
868 3 : loc = k
869 : END IF
870 : END DO
871 : ! swap symfnct
872 370 : CALL nnp_swapang(nnp%ang(i), j, loc)
873 : END DO
874 :
875 : ! sort by eta
876 : ! ang
877 370 : DO j = 1, nnp%n_ang(i) - 1
878 339 : loc = j
879 2440 : DO k = j + 1, nnp%n_ang(i)
880 2101 : IF (nnp%ang(i)%funccut(loc) == nnp%ang(i)%funccut(k) .AND. &
881 339 : nnp%ang(i)%eta(loc) > nnp%ang(i)%eta(k)) THEN
882 392 : loc = k
883 : END IF
884 : END DO
885 : ! swap symfnct
886 370 : CALL nnp_swapang(nnp%ang(i), j, loc)
887 : END DO
888 :
889 : ! sort by zeta
890 : ! ang
891 370 : DO j = 1, nnp%n_ang(i) - 1
892 339 : loc = j
893 2440 : DO k = j + 1, nnp%n_ang(i)
894 : IF (nnp%ang(i)%funccut(loc) == nnp%ang(i)%funccut(k) .AND. &
895 2101 : nnp%ang(i)%eta(loc) == nnp%ang(i)%eta(k) .AND. &
896 339 : nnp%ang(i)%zeta(loc) > nnp%ang(i)%zeta(k)) THEN
897 7 : loc = k
898 : END IF
899 : END DO
900 : ! swap symfnct
901 370 : CALL nnp_swapang(nnp%ang(i), j, loc)
902 : END DO
903 :
904 : ! sort by lambda
905 : ! ang
906 370 : DO j = 1, nnp%n_ang(i) - 1
907 339 : loc = j
908 2440 : DO k = j + 1, nnp%n_ang(i)
909 : IF (nnp%ang(i)%funccut(loc) == nnp%ang(i)%funccut(k) .AND. &
910 : nnp%ang(i)%eta(loc) == nnp%ang(i)%eta(k) .AND. &
911 2101 : nnp%ang(i)%zeta(loc) == nnp%ang(i)%zeta(k) .AND. &
912 339 : nnp%ang(i)%lam(loc) > nnp%ang(i)%lam(k)) THEN
913 148 : loc = k
914 : END IF
915 : END DO
916 : ! swap symfnct
917 370 : CALL nnp_swapang(nnp%ang(i), j, loc)
918 : END DO
919 :
920 : ! sort by ele
921 : ! ang, ele1
922 370 : DO j = 1, nnp%n_ang(i) - 1
923 339 : loc = j
924 2440 : DO k = j + 1, nnp%n_ang(i)
925 : IF (nnp%ang(i)%funccut(loc) == nnp%ang(i)%funccut(k) .AND. &
926 : nnp%ang(i)%eta(loc) == nnp%ang(i)%eta(k) .AND. &
927 : nnp%ang(i)%zeta(loc) == nnp%ang(i)%zeta(k) .AND. &
928 2101 : nnp%ang(i)%lam(loc) == nnp%ang(i)%lam(k) .AND. &
929 339 : nnp%ang(i)%nuc_ele1(loc) > nnp%ang(i)%nuc_ele1(k)) THEN
930 42 : loc = k
931 : END IF
932 : END DO
933 : ! swap symfnct
934 370 : CALL nnp_swapang(nnp%ang(i), j, loc)
935 : END DO
936 : ! ang, ele2
937 385 : DO j = 1, nnp%n_ang(i) - 1
938 339 : loc = j
939 2440 : DO k = j + 1, nnp%n_ang(i)
940 : IF (nnp%ang(i)%funccut(loc) == nnp%ang(i)%funccut(k) .AND. &
941 : nnp%ang(i)%eta(loc) == nnp%ang(i)%eta(k) .AND. &
942 : nnp%ang(i)%zeta(loc) == nnp%ang(i)%zeta(k) .AND. &
943 : nnp%ang(i)%lam(loc) == nnp%ang(i)%lam(k) .AND. &
944 2101 : nnp%ang(i)%nuc_ele1(loc) == nnp%ang(i)%nuc_ele1(k) .AND. &
945 339 : nnp%ang(i)%nuc_ele2(loc) > nnp%ang(i)%nuc_ele2(k)) THEN
946 29 : loc = k
947 : END IF
948 : END DO
949 : ! swap symfnct
950 370 : CALL nnp_swapang(nnp%ang(i), j, loc)
951 : END DO
952 : END DO
953 :
954 15 : END SUBROUTINE nnp_sort_acsf
955 :
956 : ! **************************************************************************************************
957 : !> \brief Swap two radial symmetry functions
958 : !> \param rad ...
959 : !> \param i ...
960 : !> \param j ...
961 : !> \date 2020-10-10
962 : !> \author Christoph Schran (christoph.schran@rub.de)
963 : ! **************************************************************************************************
964 1820 : SUBROUTINE nnp_swaprad(rad, i, j)
965 : TYPE(nnp_acsf_rad_type), INTENT(INOUT) :: rad
966 : INTEGER, INTENT(IN) :: i, j
967 :
968 : CHARACTER(len=2) :: tmpc
969 : INTEGER :: tmpi
970 : REAL(KIND=dp) :: tmpr
971 :
972 1820 : tmpr = rad%funccut(i)
973 1820 : rad%funccut(i) = rad%funccut(j)
974 1820 : rad%funccut(j) = tmpr
975 :
976 1820 : tmpr = rad%eta(i)
977 1820 : rad%eta(i) = rad%eta(j)
978 1820 : rad%eta(j) = tmpr
979 :
980 1820 : tmpr = rad%rs(i)
981 1820 : rad%rs(i) = rad%rs(j)
982 1820 : rad%rs(j) = tmpr
983 :
984 1820 : tmpc = rad%ele(i)
985 1820 : rad%ele(i) = rad%ele(j)
986 1820 : rad%ele(j) = tmpc
987 :
988 1820 : tmpi = rad%nuc_ele(i)
989 1820 : rad%nuc_ele(i) = rad%nuc_ele(j)
990 1820 : rad%nuc_ele(j) = tmpi
991 :
992 1820 : END SUBROUTINE nnp_swaprad
993 :
994 : ! **************************************************************************************************
995 : !> \brief Swap two angular symmetry functions
996 : !> \param ang ...
997 : !> \param i ...
998 : !> \param j ...
999 : !> \date 2020-10-10
1000 : !> \author Christoph Schran (christoph.schran@rub.de)
1001 : ! **************************************************************************************************
1002 2034 : SUBROUTINE nnp_swapang(ang, i, j)
1003 : TYPE(nnp_acsf_ang_type), INTENT(INOUT) :: ang
1004 : INTEGER, INTENT(IN) :: i, j
1005 :
1006 : CHARACTER(len=2) :: tmpc
1007 : INTEGER :: tmpi
1008 : REAL(KIND=dp) :: tmpr
1009 :
1010 2034 : tmpr = ang%funccut(i)
1011 2034 : ang%funccut(i) = ang%funccut(j)
1012 2034 : ang%funccut(j) = tmpr
1013 :
1014 2034 : tmpr = ang%eta(i)
1015 2034 : ang%eta(i) = ang%eta(j)
1016 2034 : ang%eta(j) = tmpr
1017 :
1018 2034 : tmpr = ang%zeta(i)
1019 2034 : ang%zeta(i) = ang%zeta(j)
1020 2034 : ang%zeta(j) = tmpr
1021 :
1022 2034 : tmpr = ang%prefzeta(i)
1023 2034 : ang%prefzeta(i) = ang%prefzeta(j)
1024 2034 : ang%prefzeta(j) = tmpr
1025 :
1026 2034 : tmpr = ang%lam(i)
1027 2034 : ang%lam(i) = ang%lam(j)
1028 2034 : ang%lam(j) = tmpr
1029 :
1030 2034 : tmpc = ang%ele1(i)
1031 2034 : ang%ele1(i) = ang%ele1(j)
1032 2034 : ang%ele1(j) = tmpc
1033 :
1034 2034 : tmpi = ang%nuc_ele1(i)
1035 2034 : ang%nuc_ele1(i) = ang%nuc_ele1(j)
1036 2034 : ang%nuc_ele1(j) = tmpi
1037 :
1038 2034 : tmpc = ang%ele2(i)
1039 2034 : ang%ele2(i) = ang%ele2(j)
1040 2034 : ang%ele2(j) = tmpc
1041 :
1042 2034 : tmpi = ang%nuc_ele2(i)
1043 2034 : ang%nuc_ele2(i) = ang%nuc_ele2(j)
1044 2034 : ang%nuc_ele2(j) = tmpi
1045 :
1046 2034 : END SUBROUTINE nnp_swapang
1047 :
1048 : ! **************************************************************************************************
1049 : !> \brief Initialize symmetry function groups
1050 : !> \param nnp ...
1051 : !> \date 2020-10-10
1052 : !> \author Christoph Schran (christoph.schran@rub.de)
1053 : ! **************************************************************************************************
1054 15 : SUBROUTINE nnp_init_acsf_groups(nnp)
1055 :
1056 : TYPE(nnp_type), INTENT(INOUT) :: nnp
1057 :
1058 : INTEGER :: ang, i, j, k, rad, s
1059 : REAL(KIND=dp) :: funccut
1060 :
1061 : !find out how many symmetry functions groups are needed
1062 46 : DO i = 1, nnp%n_ele
1063 31 : nnp%rad(i)%n_symfgrp = 0
1064 31 : nnp%ang(i)%n_symfgrp = 0
1065 : !search radial symmetry functions
1066 96 : DO j = 1, nnp%n_ele
1067 65 : funccut = -1.0_dp
1068 1106 : DO s = 1, nnp%n_rad(i)
1069 1075 : IF (nnp%rad(i)%ele(s) == nnp%ele(j)) THEN
1070 486 : IF (ABS(nnp%rad(i)%funccut(s) - funccut) > 1.0e-5_dp) THEN
1071 63 : nnp%rad(i)%n_symfgrp = nnp%rad(i)%n_symfgrp + 1
1072 63 : funccut = nnp%rad(i)%funccut(s)
1073 : END IF
1074 : END IF
1075 : END DO
1076 : END DO
1077 : !search angular symmetry functions
1078 111 : DO j = 1, nnp%n_ele
1079 198 : DO k = j, nnp%n_ele
1080 102 : funccut = -1.0_dp
1081 1337 : DO s = 1, nnp%n_ang(i)
1082 : IF ((nnp%ang(i)%ele1(s) == nnp%ele(j) .AND. &
1083 1170 : nnp%ang(i)%ele2(s) == nnp%ele(k)) .OR. &
1084 : (nnp%ang(i)%ele1(s) == nnp%ele(k) .AND. &
1085 102 : nnp%ang(i)%ele2(s) == nnp%ele(j))) THEN
1086 370 : IF (ABS(nnp%ang(i)%funccut(s) - funccut) > 1.0e-5_dp) THEN
1087 76 : nnp%ang(i)%n_symfgrp = nnp%ang(i)%n_symfgrp + 1
1088 76 : funccut = nnp%ang(i)%funccut(s)
1089 : END IF
1090 : END IF
1091 : END DO
1092 : END DO
1093 : END DO
1094 : END DO
1095 :
1096 : !allocate memory for symmetry functions group
1097 46 : DO i = 1, nnp%n_ele
1098 156 : ALLOCATE (nnp%rad(i)%symfgrp(nnp%rad(i)%n_symfgrp))
1099 169 : ALLOCATE (nnp%ang(i)%symfgrp(nnp%ang(i)%n_symfgrp))
1100 94 : DO j = 1, nnp%rad(i)%n_symfgrp
1101 94 : nnp%rad(i)%symfgrp(j)%n_symf = 0
1102 : END DO
1103 122 : DO j = 1, nnp%ang(i)%n_symfgrp
1104 107 : nnp%ang(i)%symfgrp(j)%n_symf = 0
1105 : END DO
1106 : END DO
1107 :
1108 : !init symmetry functions group
1109 46 : DO i = 1, nnp%n_ele
1110 : rad = 0
1111 96 : ang = 0
1112 96 : DO j = 1, nnp%n_ele
1113 65 : funccut = -1.0_dp
1114 1106 : DO s = 1, nnp%n_rad(i)
1115 1075 : IF (nnp%rad(i)%ele(s) == nnp%ele(j)) THEN
1116 486 : IF (ABS(nnp%rad(i)%funccut(s) - funccut) > 1.0e-5_dp) THEN
1117 63 : rad = rad + 1
1118 63 : funccut = nnp%rad(i)%funccut(s)
1119 63 : nnp%rad(i)%symfgrp(rad)%cutoff = funccut
1120 63 : ALLOCATE (nnp%rad(i)%symfgrp(rad)%ele(1))
1121 63 : ALLOCATE (nnp%rad(i)%symfgrp(rad)%ele_ind(1))
1122 63 : nnp%rad(i)%symfgrp(rad)%ele(1) = nnp%ele(j)
1123 63 : nnp%rad(i)%symfgrp(rad)%ele_ind(1) = j
1124 : END IF
1125 486 : nnp%rad(i)%symfgrp(rad)%n_symf = nnp%rad(i)%symfgrp(rad)%n_symf + 1
1126 : END IF
1127 : END DO
1128 : END DO
1129 111 : DO j = 1, nnp%n_ele
1130 198 : DO k = j, nnp%n_ele
1131 102 : funccut = -1.0_dp
1132 1337 : DO s = 1, nnp%n_ang(i)
1133 : IF ((nnp%ang(i)%ele1(s) == nnp%ele(j) .AND. &
1134 1170 : nnp%ang(i)%ele2(s) == nnp%ele(k)) .OR. &
1135 : (nnp%ang(i)%ele1(s) == nnp%ele(k) .AND. &
1136 102 : nnp%ang(i)%ele2(s) == nnp%ele(j))) THEN
1137 370 : IF (ABS(nnp%ang(i)%funccut(s) - funccut) > 1.0e-5_dp) THEN
1138 76 : ang = ang + 1
1139 76 : funccut = nnp%ang(i)%funccut(s)
1140 76 : nnp%ang(i)%symfgrp(ang)%cutoff = funccut
1141 76 : ALLOCATE (nnp%ang(i)%symfgrp(ang)%ele(2))
1142 76 : ALLOCATE (nnp%ang(i)%symfgrp(ang)%ele_ind(2))
1143 76 : nnp%ang(i)%symfgrp(ang)%ele(1) = nnp%ele(j)
1144 76 : nnp%ang(i)%symfgrp(ang)%ele(2) = nnp%ele(k)
1145 76 : nnp%ang(i)%symfgrp(ang)%ele_ind(1) = j
1146 76 : nnp%ang(i)%symfgrp(ang)%ele_ind(2) = k
1147 : END IF
1148 370 : nnp%ang(i)%symfgrp(ang)%n_symf = nnp%ang(i)%symfgrp(ang)%n_symf + 1
1149 : END IF
1150 : END DO
1151 : END DO
1152 : END DO
1153 : END DO
1154 :
1155 : !add symmetry functions to associated groups
1156 46 : DO i = 1, nnp%n_ele
1157 94 : DO j = 1, nnp%rad(i)%n_symfgrp
1158 189 : ALLOCATE (nnp%rad(i)%symfgrp(j)%symf(nnp%rad(i)%symfgrp(j)%n_symf))
1159 63 : rad = 0
1160 1083 : DO s = 1, nnp%n_rad(i)
1161 1052 : IF (nnp%rad(i)%ele(s) == nnp%rad(i)%symfgrp(j)%ele(1)) THEN
1162 486 : IF (ABS(nnp%rad(i)%funccut(s) - nnp%rad(i)%symfgrp(j)%cutoff) <= 1.0e-5_dp) THEN
1163 486 : rad = rad + 1
1164 486 : nnp%rad(i)%symfgrp(j)%symf(rad) = s
1165 : END IF
1166 : END IF
1167 : END DO
1168 : END DO
1169 122 : DO j = 1, nnp%ang(i)%n_symfgrp
1170 228 : ALLOCATE (nnp%ang(i)%symfgrp(j)%symf(nnp%ang(i)%symfgrp(j)%n_symf))
1171 76 : ang = 0
1172 1043 : DO s = 1, nnp%n_ang(i)
1173 : IF ((nnp%ang(i)%ele1(s) == nnp%ang(i)%symfgrp(j)%ele(1) .AND. &
1174 936 : nnp%ang(i)%ele2(s) == nnp%ang(i)%symfgrp(j)%ele(2)) .OR. &
1175 : (nnp%ang(i)%ele1(s) == nnp%ang(i)%symfgrp(j)%ele(2) .AND. &
1176 76 : nnp%ang(i)%ele2(s) == nnp%ang(i)%symfgrp(j)%ele(1))) THEN
1177 370 : IF (ABS(nnp%ang(i)%funccut(s) - nnp%ang(i)%symfgrp(j)%cutoff) <= 1.0e-5_dp) THEN
1178 370 : ang = ang + 1
1179 370 : nnp%ang(i)%symfgrp(j)%symf(ang) = s
1180 : END IF
1181 : END IF
1182 : END DO
1183 : END DO
1184 : END DO
1185 :
1186 15 : END SUBROUTINE nnp_init_acsf_groups
1187 :
1188 : ! **************************************************************************************************
1189 : !> \brief Write symmetry function information
1190 : !> \param nnp ...
1191 : !> \param para_env ...
1192 : !> \param printtag ...
1193 : !> \date 2020-10-10
1194 : !> \author Christoph Schran (christoph.schran@rub.de)
1195 : ! **************************************************************************************************
1196 15 : SUBROUTINE nnp_write_acsf(nnp, para_env, printtag)
1197 : TYPE(nnp_type), INTENT(INOUT) :: nnp
1198 : TYPE(mp_para_env_type), POINTER :: para_env
1199 : CHARACTER(LEN=*), INTENT(IN) :: printtag
1200 :
1201 : CHARACTER(len=default_string_length) :: my_label
1202 : INTEGER :: i, j, unit_nr
1203 : TYPE(cp_logger_type), POINTER :: logger
1204 :
1205 15 : NULLIFY (logger)
1206 15 : logger => cp_get_default_logger()
1207 :
1208 15 : my_label = TRIM(printtag)//"| "
1209 15 : IF (para_env%is_source()) THEN
1210 8 : unit_nr = cp_logger_get_default_unit_nr(logger)
1211 8 : WRITE (unit_nr, '(1X,A,1X,10(I2,1X))') TRIM(my_label)//" Activation functions:", nnp%actfnct(:)
1212 25 : DO i = 1, nnp%n_ele
1213 : WRITE (unit_nr, *) TRIM(my_label)//" short range atomic symmetry functions element "// &
1214 17 : nnp%ele(i)//":"
1215 279 : DO j = 1, nnp%n_rad(i)
1216 262 : 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, &
1217 262 : nnp%rad(i)%ele(j), nnp%rad(i)%eta(j), &
1218 541 : nnp%rad(i)%rs(j), nnp%rad(i)%funccut(j)
1219 : END DO
1220 220 : DO j = 1, nnp%n_ang(i)
1221 : WRITE (unit_nr, '(1X,A,1X,I3,1X,A2,1X,I2,2(1X,A2),1X,4(F6.3,1X))') &
1222 195 : TRIM(my_label), j, nnp%ele(i), 3, &
1223 195 : nnp%ang(i)%ele1(j), nnp%ang(i)%ele2(j), &
1224 195 : nnp%ang(i)%eta(j), nnp%ang(i)%lam(j), &
1225 407 : nnp%ang(i)%zeta(j), nnp%ang(i)%funccut(j)
1226 : END DO
1227 : END DO
1228 : END IF
1229 :
1230 15 : RETURN
1231 :
1232 : END SUBROUTINE nnp_write_acsf
1233 :
1234 : ! **************************************************************************************************
1235 : !> \brief Create neighbor object
1236 : !> \param nnp ...
1237 : !> \param ind ...
1238 : !> \param nat ...
1239 : !> \param neighbor ...
1240 : !> \date 2020-10-10
1241 : !> \author Christoph Schran (christoph.schran@rub.de)
1242 : ! **************************************************************************************************
1243 248177 : SUBROUTINE nnp_neighbor_create(nnp, ind, nat, neighbor)
1244 :
1245 : TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
1246 : INTEGER, INTENT(IN) :: ind, nat
1247 : TYPE(nnp_neighbor_type), INTENT(INOUT) :: neighbor
1248 :
1249 : INTEGER :: n
1250 : TYPE(cell_type), POINTER :: cell
1251 :
1252 248177 : NULLIFY (cell)
1253 248177 : CALL nnp_env_get(nnp_env=nnp, cell=cell)
1254 :
1255 248177 : CALL nnp_compute_pbc_copies(neighbor%pbc_copies, cell, nnp%max_cut)
1256 992708 : n = (SUM(neighbor%pbc_copies) + 1)*nat
1257 992708 : ALLOCATE (neighbor%dist_rad(4, n, nnp%rad(ind)%n_symfgrp))
1258 992708 : ALLOCATE (neighbor%dist_ang1(4, n, nnp%ang(ind)%n_symfgrp))
1259 992708 : ALLOCATE (neighbor%dist_ang2(4, n, nnp%ang(ind)%n_symfgrp))
1260 992708 : ALLOCATE (neighbor%ind_rad(n, nnp%rad(ind)%n_symfgrp))
1261 992708 : ALLOCATE (neighbor%ind_ang1(n, nnp%ang(ind)%n_symfgrp))
1262 992708 : ALLOCATE (neighbor%ind_ang2(n, nnp%ang(ind)%n_symfgrp))
1263 744531 : ALLOCATE (neighbor%n_rad(nnp%rad(ind)%n_symfgrp))
1264 744531 : ALLOCATE (neighbor%n_ang1(nnp%ang(ind)%n_symfgrp))
1265 744531 : ALLOCATE (neighbor%n_ang2(nnp%ang(ind)%n_symfgrp))
1266 854875 : neighbor%n_rad(:) = 0
1267 753694 : neighbor%n_ang1(:) = 0
1268 753694 : neighbor%n_ang2(:) = 0
1269 :
1270 248177 : END SUBROUTINE nnp_neighbor_create
1271 :
1272 : ! **************************************************************************************************
1273 : !> \brief Destroy neighbor object
1274 : !> \param neighbor ...
1275 : !> \date 2020-10-10
1276 : !> \author Christoph Schran (christoph.schran@rub.de)
1277 : ! **************************************************************************************************
1278 248177 : SUBROUTINE nnp_neighbor_release(neighbor)
1279 :
1280 : TYPE(nnp_neighbor_type), INTENT(INOUT) :: neighbor
1281 :
1282 248177 : DEALLOCATE (neighbor%dist_rad)
1283 248177 : DEALLOCATE (neighbor%dist_ang1)
1284 248177 : DEALLOCATE (neighbor%dist_ang2)
1285 248177 : DEALLOCATE (neighbor%ind_rad)
1286 248177 : DEALLOCATE (neighbor%ind_ang1)
1287 248177 : DEALLOCATE (neighbor%ind_ang2)
1288 248177 : DEALLOCATE (neighbor%n_rad)
1289 248177 : DEALLOCATE (neighbor%n_ang1)
1290 248177 : DEALLOCATE (neighbor%n_ang2)
1291 :
1292 248177 : END SUBROUTINE nnp_neighbor_release
1293 :
1294 : ! **************************************************************************************************
1295 : !> \brief Generate neighboring list for an atomic configuration
1296 : !> \param nnp ...
1297 : !> \param neighbor ...
1298 : !> \param i ...
1299 : !> \date 2020-10-10
1300 : !> \author Christoph Schran (christoph.schran@rub.de)
1301 : ! **************************************************************************************************
1302 248177 : SUBROUTINE nnp_compute_neighbors(nnp, neighbor, i)
1303 :
1304 : TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
1305 : TYPE(nnp_neighbor_type), INTENT(INOUT) :: neighbor
1306 : INTEGER, INTENT(IN) :: i
1307 :
1308 : INTEGER :: c1, c2, c3, ind, j, s
1309 : INTEGER, DIMENSION(3) :: nl
1310 : REAL(KIND=dp) :: norm
1311 : REAL(KIND=dp), DIMENSION(3) :: dr
1312 : TYPE(cell_type), POINTER :: cell
1313 :
1314 248177 : NULLIFY (cell)
1315 248177 : CALL nnp_env_get(nnp_env=nnp, cell=cell)
1316 :
1317 248177 : ind = nnp%ele_ind(i)
1318 :
1319 6402580 : DO j = 1, nnp%num_atoms
1320 23100087 : DO c1 = 1, 2*neighbor%pbc_copies(1) + 1
1321 16697507 : nl(1) = -neighbor%pbc_copies(1) + c1 - 1
1322 71178729 : DO c2 = 1, 2*neighbor%pbc_copies(2) + 1
1323 48326819 : nl(2) = -neighbor%pbc_copies(2) + c2 - 1
1324 208239081 : DO c3 = 1, 2*neighbor%pbc_copies(3) + 1
1325 143214755 : nl(3) = -neighbor%pbc_copies(3) + c3 - 1
1326 143214755 : IF (j == i .AND. nl(1) == 0 .AND. nl(2) == 0 .AND. nl(3) == 0) CYCLE
1327 571866312 : dr(:) = nnp%coord(:, i) - nnp%coord(:, j)
1328 : !Apply pbc, but subtract nl boxes from periodic image
1329 571866312 : dr = pbc(dr, cell, nl)
1330 571866312 : norm = NORM2(dr(:))
1331 429230766 : DO s = 1, nnp%rad(ind)%n_symfgrp
1332 429230766 : IF (nnp%ele_ind(j) == nnp%rad(ind)%symfgrp(s)%ele_ind(1)) THEN
1333 142966578 : IF (norm < nnp%rad(ind)%symfgrp(s)%cutoff) THEN
1334 3290284 : neighbor%n_rad(s) = neighbor%n_rad(s) + 1
1335 3290284 : neighbor%ind_rad(neighbor%n_rad(s), s) = j
1336 13161136 : neighbor%dist_rad(1:3, neighbor%n_rad(s), s) = dr(:)
1337 3290284 : neighbor%dist_rad(4, neighbor%n_rad(s), s) = norm
1338 : END IF
1339 : END IF
1340 : END DO
1341 524661391 : DO s = 1, nnp%ang(ind)%n_symfgrp
1342 476582749 : IF (norm < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
1343 7532482 : IF (nnp%ele_ind(j) == nnp%ang(ind)%symfgrp(s)%ele_ind(1)) THEN
1344 4042123 : neighbor%n_ang1(s) = neighbor%n_ang1(s) + 1
1345 4042123 : neighbor%ind_ang1(neighbor%n_ang1(s), s) = j
1346 16168492 : neighbor%dist_ang1(1:3, neighbor%n_ang1(s), s) = dr(:)
1347 4042123 : neighbor%dist_ang1(4, neighbor%n_ang1(s), s) = norm
1348 : END IF
1349 7532482 : IF (nnp%ele_ind(j) == nnp%ang(ind)%symfgrp(s)%ele_ind(2)) THEN
1350 2855465 : neighbor%n_ang2(s) = neighbor%n_ang2(s) + 1
1351 2855465 : neighbor%ind_ang2(neighbor%n_ang2(s), s) = j
1352 11421860 : neighbor%dist_ang2(1:3, neighbor%n_ang2(s), s) = dr(:)
1353 2855465 : neighbor%dist_ang2(4, neighbor%n_ang2(s), s) = norm
1354 : END IF
1355 : END IF
1356 : END DO
1357 : END DO
1358 : END DO
1359 : END DO
1360 : END DO
1361 :
1362 248177 : END SUBROUTINE nnp_compute_neighbors
1363 :
1364 : ! **************************************************************************************************
1365 : !> \brief Determine required pbc copies for small cells
1366 : !> \param pbc_copies ...
1367 : !> \param cell ...
1368 : !> \param cutoff ...
1369 : !> \date 2020-10-10
1370 : !> \author Christoph Schran (christoph.schran@rub.de)
1371 : ! **************************************************************************************************
1372 248177 : SUBROUTINE nnp_compute_pbc_copies(pbc_copies, cell, cutoff)
1373 : INTEGER, DIMENSION(3), INTENT(INOUT) :: pbc_copies
1374 : TYPE(cell_type), INTENT(IN), POINTER :: cell
1375 : REAL(KIND=dp), INTENT(IN) :: cutoff
1376 :
1377 : REAL(KIND=dp) :: proja, projb, projc
1378 : REAL(KIND=dp), DIMENSION(3) :: axb, axc, bxc
1379 :
1380 248177 : axb(1) = cell%hmat(2, 1)*cell%hmat(3, 2) - cell%hmat(3, 1)*cell%hmat(2, 2)
1381 248177 : axb(2) = cell%hmat(3, 1)*cell%hmat(1, 2) - cell%hmat(1, 1)*cell%hmat(3, 2)
1382 248177 : axb(3) = cell%hmat(1, 1)*cell%hmat(2, 2) - cell%hmat(2, 1)*cell%hmat(1, 2)
1383 1737239 : axb(:) = axb(:)/NORM2(axb(:))
1384 :
1385 248177 : axc(1) = cell%hmat(2, 1)*cell%hmat(3, 3) - cell%hmat(3, 1)*cell%hmat(2, 3)
1386 248177 : axc(2) = cell%hmat(3, 1)*cell%hmat(1, 3) - cell%hmat(1, 1)*cell%hmat(3, 3)
1387 248177 : axc(3) = cell%hmat(1, 1)*cell%hmat(2, 3) - cell%hmat(2, 1)*cell%hmat(1, 3)
1388 1737239 : axc(:) = axc(:)/NORM2(axc(:))
1389 :
1390 248177 : bxc(1) = cell%hmat(2, 2)*cell%hmat(3, 3) - cell%hmat(3, 2)*cell%hmat(2, 3)
1391 248177 : bxc(2) = cell%hmat(3, 2)*cell%hmat(1, 3) - cell%hmat(1, 2)*cell%hmat(3, 3)
1392 248177 : bxc(3) = cell%hmat(1, 2)*cell%hmat(2, 3) - cell%hmat(2, 2)*cell%hmat(1, 3)
1393 1737239 : bxc(:) = bxc(:)/NORM2(bxc(:))
1394 :
1395 992708 : proja = ABS(SUM(cell%hmat(:, 1)*bxc(:)))*0.5_dp
1396 992708 : projb = ABS(SUM(cell%hmat(:, 2)*axc(:)))*0.5_dp
1397 992708 : projc = ABS(SUM(cell%hmat(:, 3)*axb(:)))*0.5_dp
1398 :
1399 248177 : pbc_copies(:) = 0
1400 937730 : DO WHILE ((pbc_copies(1) + 1)*proja <= cutoff)
1401 689553 : pbc_copies(1) = pbc_copies(1) + 1
1402 : END DO
1403 937730 : DO WHILE ((pbc_copies(2) + 1)*projb <= cutoff)
1404 689553 : pbc_copies(2) = pbc_copies(2) + 1
1405 : END DO
1406 937730 : DO WHILE ((pbc_copies(3) + 1)*projc <= cutoff)
1407 689553 : pbc_copies(3) = pbc_copies(3) + 1
1408 : END DO
1409 : ! Apply non periodic setting
1410 992708 : pbc_copies(:) = pbc_copies(:)*cell%perd(:)
1411 :
1412 248177 : END SUBROUTINE nnp_compute_pbc_copies
1413 :
1414 : END MODULE nnp_acsf
|
