(git:ccc2433)
nnp_acsf.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: symtmp
71  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: forcetmp
72  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: force3tmp
73  REAL(kind=dp), DIMENSION(3) :: rvect1, rvect2, rvect3
74  TYPE(nnp_neighbor_type) :: neighbor
75 
76  CALL timeset(routinen, handle)
77 
78  !determine index of atom type
79  ind = nnp%ele_ind(i)
80 
81  ! compute neighbors of atom i
82  CALL nnp_neighbor_create(nnp, ind, nnp%num_atoms, neighbor)
83  CALL nnp_compute_neighbors(nnp, neighbor, i)
84 
85  ! Reset y:
86  nnp%rad(ind)%y = 0.0_dp
87  nnp%ang(ind)%y = 0.0_dp
88 
89  !calc forces
90  IF (PRESENT(dsymdxyz)) THEN
91  !loop over radial sym fnct grps
92  CALL timeset('nnp_acsf_radial', handle_sf)
93  DO s = 1, nnp%rad(ind)%n_symfgrp
94  ALLOCATE (symtmp(nnp%rad(ind)%symfgrp(s)%n_symf))
95  ALLOCATE (forcetmp(3, nnp%rad(ind)%symfgrp(s)%n_symf))
96  !loop over associated neighbors
97  DO j = 1, neighbor%n_rad(s)
98  rvect1 = neighbor%dist_rad(1:3, j, s)
99  r1 = neighbor%dist_rad(4, j, s)
100  CALL nnp_calc_rad(nnp, ind, s, rvect1, r1, symtmp, forcetmp)
101  jj = neighbor%ind_rad(j, s)
102  DO sf = 1, nnp%rad(ind)%symfgrp(s)%n_symf
103  m = nnp%rad(ind)%symfgrp(s)%symf(sf)
104  ! update forces into dsymdxyz
105  DO l = 1, 3
106  dsymdxyz(l, m, i) = dsymdxyz(l, m, i) + forcetmp(l, sf)
107  dsymdxyz(l, m, jj) = dsymdxyz(l, m, jj) - forcetmp(l, sf)
108  END DO
109  IF (PRESENT(stress)) THEN
110  DO l = 1, 3
111  stress(:, l, m) = stress(:, l, m) + rvect1(:)*forcetmp(l, sf)
112  END DO
113  END IF
114  nnp%rad(ind)%y(m) = nnp%rad(ind)%y(m) + symtmp(sf)
115  END DO
116  END DO
117  DEALLOCATE (symtmp)
118  DEALLOCATE (forcetmp)
119  END DO
120  CALL timestop(handle_sf)
121  !loop over angular sym fnct grps
122  CALL timeset('nnp_acsf_angular', handle_sf)
123  off = nnp%n_rad(ind)
124  DO s = 1, nnp%ang(ind)%n_symfgrp
125  ALLOCATE (symtmp(nnp%ang(ind)%symfgrp(s)%n_symf))
126  ALLOCATE (force3tmp(3, 3, nnp%ang(ind)%symfgrp(s)%n_symf))
127  !loop over associated neighbors
128  IF (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2)) THEN
129  DO j = 1, neighbor%n_ang1(s)
130  rvect1 = neighbor%dist_ang1(1:3, j, s)
131  r1 = neighbor%dist_ang1(4, j, s)
132  DO k = j + 1, neighbor%n_ang1(s)
133  rvect2 = neighbor%dist_ang1(1:3, k, s)
134  r2 = neighbor%dist_ang1(4, k, s)
135  rvect3 = rvect2 - rvect1
136  r3 = norm2(rvect3(:))
137  IF (r3 < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
138  jj = neighbor%ind_ang1(j, s)
139  kk = neighbor%ind_ang1(k, s)
140  CALL nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, &
141  r1, r2, r3, symtmp, force3tmp)
142  DO sf = 1, nnp%ang(ind)%symfgrp(s)%n_symf
143  m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
144  ! update forces into dsymdxy
145  DO l = 1, 3
146  dsymdxyz(l, m, i) = dsymdxyz(l, m, i) &
147  + force3tmp(l, 1, sf)
148  dsymdxyz(l, m, jj) = dsymdxyz(l, m, jj) &
149  + force3tmp(l, 2, sf)
150  dsymdxyz(l, m, kk) = dsymdxyz(l, m, kk) &
151  + force3tmp(l, 3, sf)
152  END DO
153  IF (PRESENT(stress)) THEN
154  DO l = 1, 3
155  stress(:, l, m) = stress(:, l, m) - rvect1(:)*force3tmp(l, 2, sf)
156  stress(:, l, m) = stress(:, l, m) - rvect2(:)*force3tmp(l, 3, sf)
157  END DO
158  END IF
159  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  DO j = 1, neighbor%n_ang1(s)
166  rvect1 = neighbor%dist_ang1(1:3, j, s)
167  r1 = neighbor%dist_ang1(4, j, s)
168  DO k = 1, neighbor%n_ang2(s)
169  rvect2 = neighbor%dist_ang2(1:3, k, s)
170  r2 = neighbor%dist_ang2(4, k, s)
171  rvect3 = rvect2 - rvect1
172  r3 = norm2(rvect3(:))
173  IF (r3 < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
174  jj = neighbor%ind_ang1(j, s)
175  kk = neighbor%ind_ang1(k, s)
176  CALL nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, &
177  r1, r2, r3, symtmp, force3tmp)
178  !loop over associated sym fncts
179  DO sf = 1, nnp%ang(ind)%symfgrp(s)%n_symf
180  m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
181  jj = neighbor%ind_ang1(j, s)
182  kk = neighbor%ind_ang2(k, s)
183  ! update forces into dsymdxy
184  DO l = 1, 3
185  dsymdxyz(l, m, i) = dsymdxyz(l, m, i) &
186  + force3tmp(l, 1, sf)
187  dsymdxyz(l, m, jj) = dsymdxyz(l, m, jj) &
188  + force3tmp(l, 2, sf)
189  dsymdxyz(l, m, kk) = dsymdxyz(l, m, kk) &
190  + force3tmp(l, 3, sf)
191  END DO
192  IF (PRESENT(stress)) THEN
193  DO l = 1, 3
194  stress(:, l, m) = stress(:, l, m) - rvect1(:)*force3tmp(l, 2, sf)
195  stress(:, l, m) = stress(:, l, m) - rvect2(:)*force3tmp(l, 3, sf)
196  END DO
197  END IF
198  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  DEALLOCATE (symtmp)
205  DEALLOCATE (force3tmp)
206  END DO
207  CALL timestop(handle_sf)
208  !no forces
209  ELSE
210  !loop over radial sym fnct grps
211  CALL timeset('nnp_acsf_radial', handle_sf)
212  DO s = 1, nnp%rad(ind)%n_symfgrp
213  ALLOCATE (symtmp(nnp%rad(ind)%symfgrp(s)%n_symf))
214  !loop over associated neighbors
215  DO j = 1, neighbor%n_rad(s)
216  rvect1 = neighbor%dist_rad(1:3, j, s)
217  r1 = neighbor%dist_rad(4, j, s)
218  CALL nnp_calc_rad(nnp, ind, s, rvect1, r1, symtmp)
219  DO sf = 1, nnp%rad(ind)%symfgrp(s)%n_symf
220  m = nnp%rad(ind)%symfgrp(s)%symf(sf)
221  nnp%rad(ind)%y(m) = nnp%rad(ind)%y(m) + symtmp(sf)
222  END DO
223  END DO
224  DEALLOCATE (symtmp)
225  END DO
226  CALL timestop(handle_sf)
227  !loop over angular sym fnct grps
228  CALL timeset('nnp_acsf_angular', handle_sf)
229  off = nnp%n_rad(ind)
230  DO s = 1, nnp%ang(ind)%n_symfgrp
231  ALLOCATE (symtmp(nnp%ang(ind)%symfgrp(s)%n_symf))
232  !loop over associated neighbors
233  IF (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2)) THEN
234  DO j = 1, neighbor%n_ang1(s)
235  rvect1 = neighbor%dist_ang1(1:3, j, s)
236  r1 = neighbor%dist_ang1(4, j, s)
237  DO k = j + 1, neighbor%n_ang1(s)
238  rvect2 = neighbor%dist_ang1(1:3, k, s)
239  r2 = neighbor%dist_ang1(4, k, s)
240  rvect3 = rvect2 - rvect1
241  r3 = norm2(rvect3(:))
242  IF (r3 < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
243  CALL nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, r1, r2, r3, symtmp)
244  DO sf = 1, nnp%ang(ind)%symfgrp(s)%n_symf
245  m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
246  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  DO j = 1, neighbor%n_ang1(s)
253  rvect1 = neighbor%dist_ang1(1:3, j, s)
254  r1 = neighbor%dist_ang1(4, j, s)
255  DO k = 1, neighbor%n_ang2(s)
256  rvect2 = neighbor%dist_ang2(1:3, k, s)
257  r2 = neighbor%dist_ang2(4, k, s)
258  rvect3 = rvect2 - rvect1
259  r3 = norm2(rvect3(:))
260  IF (r3 < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
261  CALL nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, r1, r2, r3, symtmp)
262  !loop over associated sym fncts
263  DO sf = 1, nnp%ang(ind)%symfgrp(s)%n_symf
264  m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
265  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  DEALLOCATE (symtmp)
272  END DO
273  CALL timestop(handle_sf)
274  END IF
275 
276  !check extrapolation
277  CALL nnp_check_extrapolation(nnp, ind)
278 
279  IF (PRESENT(dsymdxyz)) THEN
280  IF (PRESENT(stress)) THEN
281  CALL nnp_scale_acsf(nnp, ind, dsymdxyz, stress)
282  ELSE
283  CALL nnp_scale_acsf(nnp, ind, dsymdxyz)
284  END IF
285  ELSE
286  CALL nnp_scale_acsf(nnp, ind)
287  END IF
288 
289  CALL nnp_neighbor_release(neighbor)
290  CALL timestop(handle)
291 
292  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  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  extrapolate = nnp%output_expol
312 
313  DO j = 1, nnp%n_rad(ind)
314  IF (nnp%rad(ind)%y(j) - &
315  nnp%rad(ind)%loc_max(j) > threshold) THEN
316  extrapolate = .true.
317  ELSE IF (-nnp%rad(ind)%y(j) + &
318  nnp%rad(ind)%loc_min(j) > threshold) THEN
319  extrapolate = .true.
320  END IF
321  END DO
322  DO j = 1, nnp%n_ang(ind)
323  IF (nnp%ang(ind)%y(j) - &
324  nnp%ang(ind)%loc_max(j) > threshold) THEN
325  extrapolate = .true.
326  ELSE IF (-nnp%ang(ind)%y(j) + &
327  nnp%ang(ind)%loc_min(j) > threshold) THEN
328  extrapolate = .true.
329  END IF
330  END DO
331 
332  nnp%output_expol = extrapolate
333 
334  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  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  IF (nnp%center_acsf) THEN
355  DO j = 1, nnp%n_rad(ind)
356  nnp%arc(ind)%layer(1)%node(j) = &
357  (nnp%rad(ind)%y(j) - nnp%rad(ind)%loc_av(j))
358  END DO
359  off = nnp%n_rad(ind)
360  DO j = 1, nnp%n_ang(ind)
361  nnp%arc(ind)%layer(1)%node(j + off) = &
362  (nnp%ang(ind)%y(j) - nnp%ang(ind)%loc_av(j))
363  END DO
364 
365  IF (nnp%scale_acsf) THEN
366  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  (nnp%scmax - nnp%scmin) + nnp%scmin
371  END DO
372  off = nnp%n_rad(ind)
373  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  (nnp%scmax - nnp%scmin) + nnp%scmin
378  END DO
379  END IF
380  ELSE IF (nnp%scale_acsf) THEN
381  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  (nnp%scmax - nnp%scmin) + nnp%scmin
386  END DO
387  off = nnp%n_rad(ind)
388  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  (nnp%scmax - nnp%scmin) + nnp%scmin
393  END DO
394  ELSE IF (nnp%scale_sigma_acsf) THEN
395  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  (nnp%scmax - nnp%scmin) + nnp%scmin
400  END DO
401  off = nnp%n_rad(ind)
402  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  (nnp%scmax - nnp%scmin) + nnp%scmin
407  END DO
408  ELSE
409  DO j = 1, nnp%n_rad(ind)
410  nnp%arc(ind)%layer(1)%node(j) = nnp%rad(ind)%y(j)
411  END DO
412  off = nnp%n_rad(ind)
413  DO j = 1, nnp%n_ang(ind)
414  nnp%arc(ind)%layer(1)%node(j + off) = nnp%ang(ind)%y(j)
415  END DO
416  END IF
417 
418  IF (PRESENT(dsymdxyz)) THEN
419  IF (nnp%scale_acsf) THEN
420  DO k = 1, nnp%num_atoms
421  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  (nnp%scmax - nnp%scmin)
426  END DO
427  END DO
428  off = nnp%n_rad(ind)
429  DO k = 1, nnp%num_atoms
430  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  (nnp%scmax - nnp%scmin)
435  END DO
436  END DO
437  ELSE IF (nnp%scale_sigma_acsf) THEN
438  DO k = 1, nnp%num_atoms
439  DO j = 1, nnp%n_rad(ind)
440  dsymdxyz(:, j, k) = dsymdxyz(:, j, k)/ &
441  nnp%rad(ind)%sigma(j)* &
442  (nnp%scmax - nnp%scmin)
443  END DO
444  END DO
445  off = nnp%n_rad(ind)
446  DO k = 1, nnp%num_atoms
447  DO j = 1, nnp%n_ang(ind)
448  dsymdxyz(:, j + off, k) = dsymdxyz(:, j + off, k)/ &
449  nnp%ang(ind)%sigma(j)* &
450  (nnp%scmax - nnp%scmin)
451  END DO
452  END DO
453  END IF
454  END IF
455 
456  IF (PRESENT(stress)) THEN
457  IF (nnp%scale_acsf) THEN
458  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  (nnp%scmax - nnp%scmin)
463  END DO
464  off = nnp%n_rad(ind)
465  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  (nnp%scmax - nnp%scmin)
470  END DO
471  ELSE IF (nnp%scale_sigma_acsf) THEN
472  DO j = 1, nnp%n_rad(ind)
473  stress(:, :, j) = stress(:, :, j)/ &
474  nnp%rad(ind)%sigma(j)* &
475  (nnp%scmax - nnp%scmin)
476  END DO
477  off = nnp%n_rad(ind)
478  DO j = 1, nnp%n_ang(ind)
479  stress(:, :, j + off) = stress(:, :, j + off)/ &
480  nnp%ang(ind)%sigma(j)* &
481  (nnp%scmax - nnp%scmin)
482  END DO
483  END IF
484  END IF
485 
486  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  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  drdx = 0.0_dp
518  fcut = 0.0_dp
519  dfcutdr = 0.0_dp
520 
521  !Calculate cutoff function and partial derivative
522  funccut = nnp%rad(ind)%symfgrp(s)%cutoff !cutoff
523  SELECT CASE (nnp%cut_type)
524  CASE (nnp_cut_cos)
525  tmp = pi*r/funccut
526  fcut = 0.5_dp*(cos(tmp) + 1.0_dp)
527  IF (PRESENT(force)) THEN
528  dfcutdr = 0.5_dp*(-sin(tmp))*(pi/funccut)
529  END IF
530  CASE (nnp_cut_tanh)
531  tmp = tanh(1.0_dp - r/funccut)
532  fcut = tmp**3
533  IF (PRESENT(force)) THEN
534  dfcutdr = (-3.0_dp/funccut)*(tmp**2 - tmp**4)
535  END IF
536  CASE DEFAULT
537  cpabort("NNP| Cutoff function unknown")
538  END SELECT
539 
540  IF (PRESENT(force)) drdx(:) = rvect(:)/r
541 
542  !loop over sym fncts of sym fnct group s
543  DO sf = 1, nnp%rad(ind)%symfgrp(s)%n_symf
544  k = nnp%rad(ind)%symfgrp(s)%symf(sf) !symf indice
545  eta = nnp%rad(ind)%eta(k) !eta
546  rs = nnp%rad(ind)%rs(k) !rshift
547 
548  ! Calculate radial symmetry function
549  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  IF (PRESENT(force)) THEN
554  dsymdr = sym(sf)*(-2.0_dp*eta*(r - rs))
555  force(:, sf) = fcut*dsymdr*drdx(:) + sym(sf)*dfcutdr*drdx(:)
556  END IF
557 
558  ! combine radial symmetry function and cutoff function
559  sym(sf) = sym(sf)*fcut
560  END DO
561 
562  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  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  rsqr1 = r1**2
599  rsqr2 = r2**2
600  rsqr3 = r3**2
601 
602  !init
603  dangulardx1 = 0.0_dp
604  dangulardx2 = 0.0_dp
605  dangulardx3 = 0.0_dp
606  dr1dx = 0.0_dp
607  dr2dx = 0.0_dp
608  dr3dx = 0.0_dp
609  ftot = 0.0_dp
610  dfcutdr1 = 0.0_dp
611  dfcutdr2 = 0.0_dp
612  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  f = (rsqr3 - rsqr1 - rsqr2)
619  g = -2.0_dp*r1*r2
620  costheta = f/g
621 
622  ! Calculate cutoff function and partial derivatives
623  fcut = nnp%ang(ind)%symfgrp(s)%cutoff !cutoff
624  SELECT CASE (nnp%cut_type)
625  CASE (nnp_cut_cos)
626  tmp1 = pi*r1/fcut
627  tmp2 = pi*r2/fcut
628  tmp3 = pi*r3/fcut
629  fcut1 = 0.5_dp*(cos(tmp1) + 1.0_dp)
630  fcut2 = 0.5_dp*(cos(tmp2) + 1.0_dp)
631  fcut3 = 0.5_dp*(cos(tmp3) + 1.0_dp)
632  ftot = fcut1*fcut2*fcut3
633  IF (PRESENT(force)) THEN
634  pref = 0.5_dp*(pi/fcut)
635  dfcutdr1 = pref*(-sin(tmp1))*fcut2*fcut3
636  dfcutdr2 = pref*(-sin(tmp2))*fcut1*fcut3
637  dfcutdr3 = pref*(-sin(tmp3))*fcut1*fcut2
638  END IF
639  CASE (nnp_cut_tanh)
640  tmp1 = tanh(1.0_dp - r1/fcut)
641  tmp2 = tanh(1.0_dp - r2/fcut)
642  tmp3 = tanh(1.0_dp - r3/fcut)
643  fcut1 = tmp1**3
644  fcut2 = tmp2**3
645  fcut3 = tmp3**3
646  ftot = fcut1*fcut2*fcut3
647  IF (PRESENT(force)) THEN
648  pref = -3.0_dp/fcut
649  dfcutdr1 = pref*(tmp1**2 - tmp1**4)*fcut2*fcut3
650  dfcutdr2 = pref*(tmp2**2 - tmp2**4)*fcut1*fcut3
651  dfcutdr3 = pref*(tmp3**2 - tmp3**4)*fcut1*fcut2
652  END IF
653  CASE DEFAULT
654  cpabort("NNP| Cutoff function unknown")
655  END SELECT
656 
657  IF (PRESENT(force)) THEN
658  dr1dx(:) = rvect1(:)/r1
659  dr2dx(:) = rvect2(:)/r2
660  dr3dx(:) = rvect3(:)/r3
661  END IF
662 
663  !loop over associated sym fncts
664  DO sf = 1, nnp%ang(ind)%symfgrp(s)%n_symf
665  m = nnp%ang(ind)%symfgrp(s)%symf(sf)
666  lam = nnp%ang(ind)%lam(m) !lambda
667  zeta = nnp%ang(ind)%zeta(m) !zeta
668  prefzeta = nnp%ang(ind)%prefzeta(m) ! 2**(1-zeta)
669  eta = nnp%ang(ind)%eta(m) !eta
670 
671  tmp = (1.0_dp + lam*costheta)
672 
673  IF (tmp <= 0.0_dp) THEN
674  sym(sf) = 0.0_dp
675  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  i = nint(zeta)
683  IF (1.0_dp*i == zeta) THEN
684  tmpzeta = tmp**(i - 1) ! integer power is a LOT faster
685  ELSE
686  tmpzeta = tmp**(zeta - 1.0_dp)
687  END IF
688  angular = tmpzeta*tmp
689  ! exponential part:
690  ! exp(-eta*(r1**2+r2**2+r3**2))
691  symtmp = exp(-eta*(rsqr1 + rsqr2 + rsqr3))
692 
693  ! Partial derivatives (f/g)' = (f'g - fg')/g^2
694  IF (PRESENT(force)) THEN
695  pref = zeta*(tmpzeta)/g**2
696  DO i = 1, 3
697  dfdx1(i) = -2.0_dp*lam*(rvect1(i) + rvect2(i))
698  dfdx2(i) = 2.0_dp*lam*(rvect3(i) + rvect1(i))
699  dfdx3(i) = 2.0_dp*lam*(rvect2(i) - rvect3(i))
700 
701  tmp1 = 2.0_dp*r2*dr1dx(i)
702  tmp2 = 2.0_dp*r1*dr2dx(i)
703  dgdx1(i) = -(tmp1 + tmp2)
704  dgdx2(i) = tmp1
705  dgdx3(i) = tmp2
706 
707  dcosthetadx1(i) = dfdx1(i)*g - lam*f*dgdx1(i)
708  dcosthetadx2(i) = dfdx2(i)*g - lam*f*dgdx2(i)
709  dcosthetadx3(i) = dfdx3(i)*g - lam*f*dgdx3(i)
710 
711  dangulardx1(i) = pref*dcosthetadx1(i)
712  dangulardx2(i) = pref*dcosthetadx2(i)
713  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  pref = prefzeta
719  tmp = -2.0_dp*symtmp*eta
720  dsymdr1 = tmp*r1
721  dsymdr2 = tmp*r2
722  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  tmp = pref*symtmp*ftot
730  tmp1 = pref*angular*(ftot*dsymdr1 + dfcutdr1*symtmp)
731  tmp2 = pref*angular*(ftot*dsymdr2 + dfcutdr2*symtmp)
732  tmp3 = pref*angular*(ftot*dsymdr3 + dfcutdr3*symtmp)
733  DO i = 1, 3
734  force(i, 1, sf) = tmp*dangulardx1(i) + tmp1*dr1dx(i) + tmp2*dr2dx(i)
735  force(i, 2, sf) = tmp*dangulardx2(i) - tmp1*dr1dx(i) + tmp3*dr3dx(i)
736  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  pref = prefzeta
742  ! combine angular and exponential part with cutoff function
743  sym(sf) = pref*angular*symtmp*ftot
744  END DO
745 
746  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  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  DO i = 1, SIZE(ele)
764  CALL get_ptable_info(ele(i), number=nuc_ele(i))
765  END DO
766 
767  ! Sort both arrays
768  DO i = 1, SIZE(ele) - 1
769  minimum = nuc_ele(i)
770  loc = i
771  DO j = i + 1, SIZE(ele)
772  IF (nuc_ele(j) .LT. minimum) THEN
773  loc = j
774  minimum = nuc_ele(j)
775  END IF
776  END DO
777  tmp_nuc_ele = nuc_ele(i)
778  nuc_ele(i) = nuc_ele(loc)
779  nuc_ele(loc) = tmp_nuc_ele
780 
781  tmp_ele = ele(i)
782  ele(i) = ele(loc)
783  ele(loc) = tmp_ele
784 
785  END DO
786 
787  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  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  DO i = 1, nnp%n_ele
804  ! sort by cutoff
805  ! rad
806  DO j = 1, nnp%n_rad(i) - 1
807  loc = j
808  DO k = j + 1, nnp%n_rad(i)
809  IF (nnp%rad(i)%funccut(loc) .GT. nnp%rad(i)%funccut(k)) THEN
810  loc = k
811  END IF
812  END DO
813  ! swap symfnct
814  CALL nnp_swaprad(nnp%rad(i), j, loc)
815  END DO
816 
817  ! sort by eta
818  ! rad
819  DO j = 1, nnp%n_rad(i) - 1
820  loc = j
821  DO k = j + 1, nnp%n_rad(i)
822  IF (nnp%rad(i)%funccut(loc) .EQ. nnp%rad(i)%funccut(k) .AND. &
823  nnp%rad(i)%eta(loc) .GT. nnp%rad(i)%eta(k)) THEN
824  loc = k
825  END IF
826  END DO
827  ! swap symfnct
828  CALL nnp_swaprad(nnp%rad(i), j, loc)
829  END DO
830 
831  ! sort by rshift
832  ! rad
833  DO j = 1, nnp%n_rad(i) - 1
834  loc = j
835  DO k = j + 1, nnp%n_rad(i)
836  IF (nnp%rad(i)%funccut(loc) .EQ. nnp%rad(i)%funccut(k) .AND. &
837  nnp%rad(i)%eta(loc) .EQ. nnp%rad(i)%eta(k) .AND. &
838  nnp%rad(i)%rs(loc) .GT. nnp%rad(i)%rs(k)) THEN
839  loc = k
840  END IF
841  END DO
842  ! swap symfnct
843  CALL nnp_swaprad(nnp%rad(i), j, loc)
844  END DO
845 
846  ! sort by ele
847  ! rad
848  DO j = 1, nnp%n_rad(i) - 1
849  loc = j
850  DO k = j + 1, nnp%n_rad(i)
851  IF (nnp%rad(i)%funccut(loc) .EQ. nnp%rad(i)%funccut(k) .AND. &
852  nnp%rad(i)%eta(loc) .EQ. nnp%rad(i)%eta(k) .AND. &
853  nnp%rad(i)%rs(loc) .EQ. nnp%rad(i)%rs(k) .AND. &
854  nnp%rad(i)%nuc_ele(loc) .GT. nnp%rad(i)%nuc_ele(k)) THEN
855  loc = k
856  END IF
857  END DO
858  ! swap symfnct
859  CALL nnp_swaprad(nnp%rad(i), j, loc)
860  END DO
861 
862  ! ang
863  ! sort by cutoff
864  DO j = 1, nnp%n_ang(i) - 1
865  loc = j
866  DO k = j + 1, nnp%n_ang(i)
867  IF (nnp%ang(i)%funccut(loc) .GT. nnp%ang(i)%funccut(k)) THEN
868  loc = k
869  END IF
870  END DO
871  ! swap symfnct
872  CALL nnp_swapang(nnp%ang(i), j, loc)
873  END DO
874 
875  ! sort by eta
876  ! ang
877  DO j = 1, nnp%n_ang(i) - 1
878  loc = j
879  DO k = j + 1, nnp%n_ang(i)
880  IF (nnp%ang(i)%funccut(loc) .EQ. nnp%ang(i)%funccut(k) .AND. &
881  nnp%ang(i)%eta(loc) .GT. nnp%ang(i)%eta(k)) THEN
882  loc = k
883  END IF
884  END DO
885  ! swap symfnct
886  CALL nnp_swapang(nnp%ang(i), j, loc)
887  END DO
888 
889  ! sort by zeta
890  ! ang
891  DO j = 1, nnp%n_ang(i) - 1
892  loc = j
893  DO k = j + 1, nnp%n_ang(i)
894  IF (nnp%ang(i)%funccut(loc) .EQ. nnp%ang(i)%funccut(k) .AND. &
895  nnp%ang(i)%eta(loc) .EQ. nnp%ang(i)%eta(k) .AND. &
896  nnp%ang(i)%zeta(loc) .GT. nnp%ang(i)%zeta(k)) THEN
897  loc = k
898  END IF
899  END DO
900  ! swap symfnct
901  CALL nnp_swapang(nnp%ang(i), j, loc)
902  END DO
903 
904  ! sort by lambda
905  ! ang
906  DO j = 1, nnp%n_ang(i) - 1
907  loc = j
908  DO k = j + 1, nnp%n_ang(i)
909  IF (nnp%ang(i)%funccut(loc) .EQ. nnp%ang(i)%funccut(k) .AND. &
910  nnp%ang(i)%eta(loc) .EQ. nnp%ang(i)%eta(k) .AND. &
911  nnp%ang(i)%zeta(loc) .EQ. nnp%ang(i)%zeta(k) .AND. &
912  nnp%ang(i)%lam(loc) .GT. nnp%ang(i)%lam(k)) THEN
913  loc = k
914  END IF
915  END DO
916  ! swap symfnct
917  CALL nnp_swapang(nnp%ang(i), j, loc)
918  END DO
919 
920  ! sort by ele
921  ! ang, ele1
922  DO j = 1, nnp%n_ang(i) - 1
923  loc = j
924  DO k = j + 1, nnp%n_ang(i)
925  IF (nnp%ang(i)%funccut(loc) .EQ. nnp%ang(i)%funccut(k) .AND. &
926  nnp%ang(i)%eta(loc) .EQ. nnp%ang(i)%eta(k) .AND. &
927  nnp%ang(i)%zeta(loc) .EQ. nnp%ang(i)%zeta(k) .AND. &
928  nnp%ang(i)%lam(loc) .EQ. nnp%ang(i)%lam(k) .AND. &
929  nnp%ang(i)%nuc_ele1(loc) .GT. nnp%ang(i)%nuc_ele1(k)) THEN
930  loc = k
931  END IF
932  END DO
933  ! swap symfnct
934  CALL nnp_swapang(nnp%ang(i), j, loc)
935  END DO
936  ! ang, ele2
937  DO j = 1, nnp%n_ang(i) - 1
938  loc = j
939  DO k = j + 1, nnp%n_ang(i)
940  IF (nnp%ang(i)%funccut(loc) .EQ. nnp%ang(i)%funccut(k) .AND. &
941  nnp%ang(i)%eta(loc) .EQ. nnp%ang(i)%eta(k) .AND. &
942  nnp%ang(i)%zeta(loc) .EQ. nnp%ang(i)%zeta(k) .AND. &
943  nnp%ang(i)%lam(loc) .EQ. nnp%ang(i)%lam(k) .AND. &
944  nnp%ang(i)%nuc_ele1(loc) .EQ. nnp%ang(i)%nuc_ele1(k) .AND. &
945  nnp%ang(i)%nuc_ele2(loc) .GT. nnp%ang(i)%nuc_ele2(k)) THEN
946  loc = k
947  END IF
948  END DO
949  ! swap symfnct
950  CALL nnp_swapang(nnp%ang(i), j, loc)
951  END DO
952  END DO
953 
954  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  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  tmpr = rad%funccut(i)
973  rad%funccut(i) = rad%funccut(j)
974  rad%funccut(j) = tmpr
975 
976  tmpr = rad%eta(i)
977  rad%eta(i) = rad%eta(j)
978  rad%eta(j) = tmpr
979 
980  tmpr = rad%rs(i)
981  rad%rs(i) = rad%rs(j)
982  rad%rs(j) = tmpr
983 
984  tmpc = rad%ele(i)
985  rad%ele(i) = rad%ele(j)
986  rad%ele(j) = tmpc
987 
988  tmpi = rad%nuc_ele(i)
989  rad%nuc_ele(i) = rad%nuc_ele(j)
990  rad%nuc_ele(j) = tmpi
991 
992  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  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  tmpr = ang%funccut(i)
1011  ang%funccut(i) = ang%funccut(j)
1012  ang%funccut(j) = tmpr
1013 
1014  tmpr = ang%eta(i)
1015  ang%eta(i) = ang%eta(j)
1016  ang%eta(j) = tmpr
1017 
1018  tmpr = ang%zeta(i)
1019  ang%zeta(i) = ang%zeta(j)
1020  ang%zeta(j) = tmpr
1021 
1022  tmpr = ang%prefzeta(i)
1023  ang%prefzeta(i) = ang%prefzeta(j)
1024  ang%prefzeta(j) = tmpr
1025 
1026  tmpr = ang%lam(i)
1027  ang%lam(i) = ang%lam(j)
1028  ang%lam(j) = tmpr
1029 
1030  tmpc = ang%ele1(i)
1031  ang%ele1(i) = ang%ele1(j)
1032  ang%ele1(j) = tmpc
1033 
1034  tmpi = ang%nuc_ele1(i)
1035  ang%nuc_ele1(i) = ang%nuc_ele1(j)
1036  ang%nuc_ele1(j) = tmpi
1037 
1038  tmpc = ang%ele2(i)
1039  ang%ele2(i) = ang%ele2(j)
1040  ang%ele2(j) = tmpc
1041 
1042  tmpi = ang%nuc_ele2(i)
1043  ang%nuc_ele2(i) = ang%nuc_ele2(j)
1044  ang%nuc_ele2(j) = tmpi
1045 
1046  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  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  DO i = 1, nnp%n_ele
1063  nnp%rad(i)%n_symfgrp = 0
1064  nnp%ang(i)%n_symfgrp = 0
1065  !search radial symmetry functions
1066  DO j = 1, nnp%n_ele
1067  funccut = -1.0_dp
1068  DO s = 1, nnp%n_rad(i)
1069  IF (nnp%rad(i)%ele(s) == nnp%ele(j)) THEN
1070  IF (abs(nnp%rad(i)%funccut(s) - funccut) > 1.0e-5_dp) THEN
1071  nnp%rad(i)%n_symfgrp = nnp%rad(i)%n_symfgrp + 1
1072  funccut = nnp%rad(i)%funccut(s)
1073  END IF
1074  END IF
1075  END DO
1076  END DO
1077  !search angular symmetry functions
1078  DO j = 1, nnp%n_ele
1079  DO k = j, nnp%n_ele
1080  funccut = -1.0_dp
1081  DO s = 1, nnp%n_ang(i)
1082  IF ((nnp%ang(i)%ele1(s) == nnp%ele(j) .AND. &
1083  nnp%ang(i)%ele2(s) == nnp%ele(k)) .OR. &
1084  (nnp%ang(i)%ele1(s) == nnp%ele(k) .AND. &
1085  nnp%ang(i)%ele2(s) == nnp%ele(j))) THEN
1086  IF (abs(nnp%ang(i)%funccut(s) - funccut) > 1.0e-5_dp) THEN
1087  nnp%ang(i)%n_symfgrp = nnp%ang(i)%n_symfgrp + 1
1088  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  DO i = 1, nnp%n_ele
1098  ALLOCATE (nnp%rad(i)%symfgrp(nnp%rad(i)%n_symfgrp))
1099  ALLOCATE (nnp%ang(i)%symfgrp(nnp%ang(i)%n_symfgrp))
1100  DO j = 1, nnp%rad(i)%n_symfgrp
1101  nnp%rad(i)%symfgrp(j)%n_symf = 0
1102  END DO
1103  DO j = 1, nnp%ang(i)%n_symfgrp
1104  nnp%ang(i)%symfgrp(j)%n_symf = 0
1105  END DO
1106  END DO
1107 
1108  !init symmetry functions group
1109  DO i = 1, nnp%n_ele
1110  rad = 0
1111  ang = 0
1112  DO j = 1, nnp%n_ele
1113  funccut = -1.0_dp
1114  DO s = 1, nnp%n_rad(i)
1115  IF (nnp%rad(i)%ele(s) == nnp%ele(j)) THEN
1116  IF (abs(nnp%rad(i)%funccut(s) - funccut) > 1.0e-5_dp) THEN
1117  rad = rad + 1
1118  funccut = nnp%rad(i)%funccut(s)
1119  nnp%rad(i)%symfgrp(rad)%cutoff = funccut
1120  ALLOCATE (nnp%rad(i)%symfgrp(rad)%ele(1))
1121  ALLOCATE (nnp%rad(i)%symfgrp(rad)%ele_ind(1))
1122  nnp%rad(i)%symfgrp(rad)%ele(1) = nnp%ele(j)
1123  nnp%rad(i)%symfgrp(rad)%ele_ind(1) = j
1124  END IF
1125  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  DO j = 1, nnp%n_ele
1130  DO k = j, nnp%n_ele
1131  funccut = -1.0_dp
1132  DO s = 1, nnp%n_ang(i)
1133  IF ((nnp%ang(i)%ele1(s) == nnp%ele(j) .AND. &
1134  nnp%ang(i)%ele2(s) == nnp%ele(k)) .OR. &
1135  (nnp%ang(i)%ele1(s) == nnp%ele(k) .AND. &
1136  nnp%ang(i)%ele2(s) == nnp%ele(j))) THEN
1137  IF (abs(nnp%ang(i)%funccut(s) - funccut) > 1.0e-5_dp) THEN
1138  ang = ang + 1
1139  funccut = nnp%ang(i)%funccut(s)
1140  nnp%ang(i)%symfgrp(ang)%cutoff = funccut
1141  ALLOCATE (nnp%ang(i)%symfgrp(ang)%ele(2))
1142  ALLOCATE (nnp%ang(i)%symfgrp(ang)%ele_ind(2))
1143  nnp%ang(i)%symfgrp(ang)%ele(1) = nnp%ele(j)
1144  nnp%ang(i)%symfgrp(ang)%ele(2) = nnp%ele(k)
1145  nnp%ang(i)%symfgrp(ang)%ele_ind(1) = j
1146  nnp%ang(i)%symfgrp(ang)%ele_ind(2) = k
1147  END IF
1148  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  DO i = 1, nnp%n_ele
1157  DO j = 1, nnp%rad(i)%n_symfgrp
1158  ALLOCATE (nnp%rad(i)%symfgrp(j)%symf(nnp%rad(i)%symfgrp(j)%n_symf))
1159  rad = 0
1160  DO s = 1, nnp%n_rad(i)
1161  IF (nnp%rad(i)%ele(s) == nnp%rad(i)%symfgrp(j)%ele(1)) THEN
1162  IF (abs(nnp%rad(i)%funccut(s) - nnp%rad(i)%symfgrp(j)%cutoff) <= 1.0e-5_dp) THEN
1163  rad = rad + 1
1164  nnp%rad(i)%symfgrp(j)%symf(rad) = s
1165  END IF
1166  END IF
1167  END DO
1168  END DO
1169  DO j = 1, nnp%ang(i)%n_symfgrp
1170  ALLOCATE (nnp%ang(i)%symfgrp(j)%symf(nnp%ang(i)%symfgrp(j)%n_symf))
1171  ang = 0
1172  DO s = 1, nnp%n_ang(i)
1173  IF ((nnp%ang(i)%ele1(s) == nnp%ang(i)%symfgrp(j)%ele(1) .AND. &
1174  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  nnp%ang(i)%ele2(s) == nnp%ang(i)%symfgrp(j)%ele(1))) THEN
1177  IF (abs(nnp%ang(i)%funccut(s) - nnp%ang(i)%symfgrp(j)%cutoff) <= 1.0e-5_dp) THEN
1178  ang = ang + 1
1179  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  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  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  NULLIFY (logger)
1206  logger => cp_get_default_logger()
1207 
1208  my_label = trim(printtag)//"| "
1209  IF (para_env%is_source()) THEN
1210  unit_nr = cp_logger_get_default_unit_nr(logger)
1211  WRITE (unit_nr, '(1X,A,1X,10(I2,1X))') trim(my_label)//" Activation functions:", nnp%actfnct(:)
1212  DO i = 1, nnp%n_ele
1213  WRITE (unit_nr, *) trim(my_label)//" short range atomic symmetry functions element "// &
1214  nnp%ele(i)//":"
1215  DO j = 1, nnp%n_rad(i)
1216  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  nnp%rad(i)%ele(j), nnp%rad(i)%eta(j), &
1218  nnp%rad(i)%rs(j), nnp%rad(i)%funccut(j)
1219  END DO
1220  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  trim(my_label), j, nnp%ele(i), 3, &
1223  nnp%ang(i)%ele1(j), nnp%ang(i)%ele2(j), &
1224  nnp%ang(i)%eta(j), nnp%ang(i)%lam(j), &
1225  nnp%ang(i)%zeta(j), nnp%ang(i)%funccut(j)
1226  END DO
1227  END DO
1228  END IF
1229 
1230  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  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  NULLIFY (cell)
1253  CALL nnp_env_get(nnp_env=nnp, cell=cell)
1254 
1255  CALL nnp_compute_pbc_copies(neighbor%pbc_copies, cell, nnp%max_cut)
1256  n = (sum(neighbor%pbc_copies) + 1)*nat
1257  ALLOCATE (neighbor%dist_rad(4, n, nnp%rad(ind)%n_symfgrp))
1258  ALLOCATE (neighbor%dist_ang1(4, n, nnp%ang(ind)%n_symfgrp))
1259  ALLOCATE (neighbor%dist_ang2(4, n, nnp%ang(ind)%n_symfgrp))
1260  ALLOCATE (neighbor%ind_rad(n, nnp%rad(ind)%n_symfgrp))
1261  ALLOCATE (neighbor%ind_ang1(n, nnp%ang(ind)%n_symfgrp))
1262  ALLOCATE (neighbor%ind_ang2(n, nnp%ang(ind)%n_symfgrp))
1263  ALLOCATE (neighbor%n_rad(nnp%rad(ind)%n_symfgrp))
1264  ALLOCATE (neighbor%n_ang1(nnp%ang(ind)%n_symfgrp))
1265  ALLOCATE (neighbor%n_ang2(nnp%ang(ind)%n_symfgrp))
1266  neighbor%n_rad(:) = 0
1267  neighbor%n_ang1(:) = 0
1268  neighbor%n_ang2(:) = 0
1269 
1270  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  SUBROUTINE nnp_neighbor_release(neighbor)
1279 
1280  TYPE(nnp_neighbor_type), INTENT(INOUT) :: neighbor
1281 
1282  DEALLOCATE (neighbor%dist_rad)
1283  DEALLOCATE (neighbor%dist_ang1)
1284  DEALLOCATE (neighbor%dist_ang2)
1285  DEALLOCATE (neighbor%ind_rad)
1286  DEALLOCATE (neighbor%ind_ang1)
1287  DEALLOCATE (neighbor%ind_ang2)
1288  DEALLOCATE (neighbor%n_rad)
1289  DEALLOCATE (neighbor%n_ang1)
1290  DEALLOCATE (neighbor%n_ang2)
1291 
1292  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  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  NULLIFY (cell)
1315  CALL nnp_env_get(nnp_env=nnp, cell=cell)
1316 
1317  ind = nnp%ele_ind(i)
1318 
1319  DO j = 1, nnp%num_atoms
1320  DO c1 = 1, 2*neighbor%pbc_copies(1) + 1
1321  nl(1) = -neighbor%pbc_copies(1) + c1 - 1
1322  DO c2 = 1, 2*neighbor%pbc_copies(2) + 1
1323  nl(2) = -neighbor%pbc_copies(2) + c2 - 1
1324  DO c3 = 1, 2*neighbor%pbc_copies(3) + 1
1325  nl(3) = -neighbor%pbc_copies(3) + c3 - 1
1326  IF (j == i .AND. nl(1) == 0 .AND. nl(2) == 0 .AND. nl(3) == 0) cycle
1327  dr(:) = nnp%coord(:, i) - nnp%coord(:, j)
1328  !Apply pbc, but subtract nl boxes from periodic image
1329  dr = pbc(dr, cell, nl)
1330  norm = norm2(dr(:))
1331  DO s = 1, nnp%rad(ind)%n_symfgrp
1332  IF (nnp%ele_ind(j) == nnp%rad(ind)%symfgrp(s)%ele_ind(1)) THEN
1333  IF (norm < nnp%rad(ind)%symfgrp(s)%cutoff) THEN
1334  neighbor%n_rad(s) = neighbor%n_rad(s) + 1
1335  neighbor%ind_rad(neighbor%n_rad(s), s) = j
1336  neighbor%dist_rad(1:3, neighbor%n_rad(s), s) = dr(:)
1337  neighbor%dist_rad(4, neighbor%n_rad(s), s) = norm
1338  END IF
1339  END IF
1340  END DO
1341  DO s = 1, nnp%ang(ind)%n_symfgrp
1342  IF (norm < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
1343  IF (nnp%ele_ind(j) == nnp%ang(ind)%symfgrp(s)%ele_ind(1)) THEN
1344  neighbor%n_ang1(s) = neighbor%n_ang1(s) + 1
1345  neighbor%ind_ang1(neighbor%n_ang1(s), s) = j
1346  neighbor%dist_ang1(1:3, neighbor%n_ang1(s), s) = dr(:)
1347  neighbor%dist_ang1(4, neighbor%n_ang1(s), s) = norm
1348  END IF
1349  IF (nnp%ele_ind(j) == nnp%ang(ind)%symfgrp(s)%ele_ind(2)) THEN
1350  neighbor%n_ang2(s) = neighbor%n_ang2(s) + 1
1351  neighbor%ind_ang2(neighbor%n_ang2(s), s) = j
1352  neighbor%dist_ang2(1:3, neighbor%n_ang2(s), s) = dr(:)
1353  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  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  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  axb(1) = cell%hmat(2, 1)*cell%hmat(3, 2) - cell%hmat(3, 1)*cell%hmat(2, 2)
1381  axb(2) = cell%hmat(3, 1)*cell%hmat(1, 2) - cell%hmat(1, 1)*cell%hmat(3, 2)
1382  axb(3) = cell%hmat(1, 1)*cell%hmat(2, 2) - cell%hmat(2, 1)*cell%hmat(1, 2)
1383  axb(:) = axb(:)/norm2(axb(:))
1384 
1385  axc(1) = cell%hmat(2, 1)*cell%hmat(3, 3) - cell%hmat(3, 1)*cell%hmat(2, 3)
1386  axc(2) = cell%hmat(3, 1)*cell%hmat(1, 3) - cell%hmat(1, 1)*cell%hmat(3, 3)
1387  axc(3) = cell%hmat(1, 1)*cell%hmat(2, 3) - cell%hmat(2, 1)*cell%hmat(1, 3)
1388  axc(:) = axc(:)/norm2(axc(:))
1389 
1390  bxc(1) = cell%hmat(2, 2)*cell%hmat(3, 3) - cell%hmat(3, 2)*cell%hmat(2, 3)
1391  bxc(2) = cell%hmat(3, 2)*cell%hmat(1, 3) - cell%hmat(1, 2)*cell%hmat(3, 3)
1392  bxc(3) = cell%hmat(1, 2)*cell%hmat(2, 3) - cell%hmat(2, 2)*cell%hmat(1, 3)
1393  bxc(:) = bxc(:)/norm2(bxc(:))
1394 
1395  proja = abs(sum(cell%hmat(:, 1)*bxc(:)))*0.5_dp
1396  projb = abs(sum(cell%hmat(:, 2)*axc(:)))*0.5_dp
1397  projc = abs(sum(cell%hmat(:, 3)*axb(:)))*0.5_dp
1398 
1399  pbc_copies(:) = 0
1400  DO WHILE ((pbc_copies(1) + 1)*proja <= cutoff)
1401  pbc_copies(1) = pbc_copies(1) + 1
1402  END DO
1403  DO WHILE ((pbc_copies(2) + 1)*projb <= cutoff)
1404  pbc_copies(2) = pbc_copies(2) + 1
1405  END DO
1406  DO WHILE ((pbc_copies(3) + 1)*projc <= cutoff)
1407  pbc_copies(3) = pbc_copies(3) + 1
1408  END DO
1409  ! Apply non periodic setting
1410  pbc_copies(:) = pbc_copies(:)*cell%perd(:)
1411 
1412  END SUBROUTINE nnp_compute_pbc_copies
1413 
1414 END MODULE nnp_acsf
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition: cp_log_handling.F:567
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
nnp_acsf
Functionality for atom centered symmetry functions for neural network potentials.
Definition: nnp_acsf.F:14
nnp_acsf::nnp_calc_acsf
subroutine, public nnp_calc_acsf(nnp, i, dsymdxyz, stress)
Calculate atom centered symmetry functions for given atom i.
Definition: nnp_acsf.F:59
nnp_acsf::nnp_write_acsf
subroutine, public nnp_write_acsf(nnp, para_env, printtag)
Write symmetry function information.
Definition: nnp_acsf.F:1197
nnp_acsf::nnp_init_acsf_groups
subroutine, public nnp_init_acsf_groups(nnp)
Initialize symmetry function groups.
Definition: nnp_acsf.F:1055
nnp_acsf::nnp_sort_acsf
subroutine, public nnp_sort_acsf(nnp)
Sort symmetry functions according to different criteria.
Definition: nnp_acsf.F:796
nnp_acsf::nnp_sort_ele
subroutine, public nnp_sort_ele(ele, nuc_ele)
Sort element array according to atomic number.
Definition: nnp_acsf.F:756
nnp_environment_types
Data types for neural network potentials.
Definition: nnp_environment_types.F:13
nnp_environment_types::nnp_cut_tanh
integer, parameter, public nnp_cut_tanh
Definition: nnp_environment_types.F:63
nnp_environment_types::nnp_cut_cos
integer, parameter, public nnp_cut_cos
derived data types
Definition: nnp_environment_types.F:63
nnp_environment_types::nnp_env_get
subroutine, public nnp_env_get(nnp_env, nnp_forces, subsys, atomic_kind_set, particle_set, local_particles, molecule_kind_set, molecule_set, local_molecules, nnp_input, force_env_input, cell, cell_ref, use_ref_cell, nnp_potential_energy, virial)
Returns various attributes of the nnp environment.
Definition: nnp_environment_types.F:418
periodic_table
Periodic Table related data definitions.
Definition: periodic_table.F:28
periodic_table::get_ptable_info
subroutine, public get_ptable_info(symbol, number, amass, ielement, covalent_radius, metallic_radius, vdw_radius, found)
Pass information about the kind given the element symbol.
Definition: periodic_table.F:85